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STATEMENT  OF  PROBLEM 


The  general  problem  being  addressed  in  this  investigation  is  how  to  develop  accurate  theoretical 
models  for  microwave  and  millimeter- wave  propagation  and  scattering  in  vegetation.  Such  models  are  needed 
since  vegetation  along  the  propagation  channel  contributes  to  many  problems  and  limitations  that  affect 
mobile  radio  systems.  Of  particular  interest  to  the  US  Army  is  the  development  of  models  that  provide  real 
time  data  to  the  soldier  in  the  field.  Hence,  both  rigorous  and  approximate  models  are  of  interest. 

At  frequencies  above  3GHz,  the  wavelength  is  of  the  order  of,  or  small  compared  to,  the  dimensions 
and  spacings  of  vegetation,  such  as  leaves,  twigs  and  branches.  Vegetation  has  to  be  modeled  as  a  random 
medium  of  discrete,  lossy  scatterers.  Multi  scattering  is  important  and  both  the  coherent  and  incoherent  field 
components  have  to  be  taken  into  account.  To  account  for  these  characteristics,  the  theoretical  models 
developed  here  for  propagation  and  scattering  in  vegetation  at  high  frequency  are  based  on  radiative  energy 
transfer  theory  (transport  theory). 

The  specific  problem  analyzed  is  that  of  a  periodic  sequence  of  Gaussian  pulses  incident  from  free 
space  (air)  onto  a  forest  region  (vegetation).  The  forest  is  modeled  as  a  half-space  of  randomly  distributed 
particles  that  scatter  and  absorb  electromagnetic  energy.  The  incident  pulse  train  is  taken  to  be  either  a 
collimated  (cylindrical)  wave  beam,  which  is  pertinent  for  IR  and  optical  communications,  or  a  spherical  wave 
that  is  confined  within  a  certain  solid  angle,  which  is  characteristic  of  radiation  produced  by  a  microwave  or 


mm- wave  antenna. 


SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS 


The  major  contribution  of  this  study  is  the  development  of  several  new  theories  -  both  rigorous  and 
approximate  -  capable  of  characterizing  the  propagation  and  scattering  of  microwave  and  millimeter-wave 
pulsed  beam  waves  in  vegetation.  These  theories  allowed  for  a  comprehensive  characterization  of  the 
influence  of  vegetation  on  the  propagation  of  beam  waves,  which  included  a  description  of  their  attenuation, 
their  angular  spread,  their  distortion  due  to  pulse  broadening,  and  the  determination  of  out-of-the-beam 
scattering. 

The  theories  are  based  on  solutions  to  the  equation  of  radiative  transfer  in  cylindrical  coordinates. 
Such  theories  are  new  and  provide  better  understanding  of  the  effects  of  vegetation  in  the  propagation 
channel.  Two  cases  are  considered.  The  first  case  deals  with  an  incident  collimated  beam  wave  pulse  train 
that  enters  a  forest  half- space.  The  second  case  assumes  that  the  incident  radiation  upon  the  forest  is  a 
diverging  spherical  beam  wave.  In  both  cases,  vegetation  is  modeled  as  a  random  medium  of  discrete 
scatterers. 

To  allow  for  real  time  determination  of  the  attenuation  and  scattering  of  the  field  in  vegetation,  three 
approximate  theories  for  both  the  collimated  and  spherical  beams  are  developed.  The  first  theory,  a  first  order 
multiple  scattering  theory,  allows  estimation  of  the  intensity  at  small  penetration  depths  without  resorting  to 
computer  intensive  calculations  while  maintaining  good  accuracy.  The  second  new  theory  yields  an 
asymptotic  solution  to  the  transport  equation,  which  is  valid  at  large  penetration  depths  into  the  forest.  This 
allows  a  characterization  of  the  field  deep  in  the  forest  also  without  extensive  calculations.  The  third  theory, 
which  combines  the  latter  two  approximate  theories,  attempts  to  close  the  gap  between  the  first  order  theory 
and  the  asymptotic  theory.  As  opposed  to  the  rigorous  theories,  these  approximate  theories  have  been  found  to 
be  numerically  efficient.  The  importance  of  these  approximate  theories  is  that  they  provide  real-time 
information  to  the  soldier  in  the  field  quickly  as  opposed  to  the  rigorous  theory  which,  being  computationally 
intensive,  is  time  consuming. 
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1.  INTRODUCTION 


For  line-of-sight  communication,  cellular  communication  in  particular,  current  interest  centers  on 
radio-link  performance,  and  how  it  is  affected  by  wave  attenuation,  fading  and  co-channel  interference.  When 
vegetation,  such  as  a  forest,  lies  along  the  path  of  a  radio-link,  the  radio  performance  is  affected  by  strong 
multiscattering  effects. 

There  are  two  methods  that  are  usually  used  to  study  multiscattering  effects  in  random  media1, 
namely,  analytical  theory  and  transport  theory  [1].  Analytical  theory  is  a  very  rigorous  mathematical 
approach  based  on  Maxwell’s  Equations.  It  is  very  complex  and  obtaining  solutions  often  requires 
introducing  strong  simplifications  which  limit  the  solutions  obtained  to  restricted  parameter  ranges.  In 
contrast,  radiative  transfer  theory  deals  with  the  transfer  of  energy  through  the  multiscattering  medium.  In  this 
theory,  the  basic  equation  that  is  solved  is  the  equation  of  radiative  transfer  or  transport  equation.  The 
radiative  transfer  theory  developed  heuristically  from  the  conservation  of  energy  principle  in  radiation  space. 
The  transport  equation  is  equivalent  to  Boltzmann’s  equation  found  in  the  kinetic  theory  of  gases  and  in 
neutron  transport  theory  and  is  less  rigorous  than  the  analytical  theory.  However,  transport  theory  has  been 
very  successfully  applied  in  the  study  of  many  radiation  problems,  such  as  optical  propagation  through  the 
atmosphere,  remote  sensing  and  radiation  from  stars. 

In  previous  work,  continuous  wave  (CW)  millimeter  wave  and  plane  wave  pulse  propagation  in 
vegetation  were  studied  using  the  scalar  transport  theory  [2-6].  In  these  studies,  interest  focused  on  the 
determination  of  the  range  and  directional  dependency  of  the  received  power  as  well  as  on  pulse  broadening 
and  distortion.  The  scalar  transport  equation  is  capable  of  specifying  the  total  energy  density  of  radiation  in 
two  orthogonal  polarizations,  but  not  polarization  or  depolarization  effects  (see  [4]  for  experimental 
justification  of  their  neglect  in  these  studies  for  millimeter  waves  in  vegetation).  In  the  earlier  developed 
theory  of  a  plane  wave  incident  upon  the  forest  half- space,  it  was  shown  that  the  range  dependence  in  the 
forest  (treated  as  a  random  medium)  is  not  simply  an  exponential  decrease  at  constant  attenuation  rate.  What 
actually  occurs  for  the  received  power  is  a  high  attenuation  rate  at  short  distances  into  the  medium  that 
evolves  into  a  much  lower  attenuation  rate  at  large  distances.  The  theory  explains  this  in  terms  of  the 
interaction  between  the  coherent  and  incoherent  field  components.  The  coherent  component,  dominating  at 

1  In  this  study,  a  forest  half-space  is  modeled  as  a  random  medium. 
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short  distances,  is  highly  attenuated  (by  absorption  and  scattering  )  while  the  incoherent  component,  which  is 
generated  by  the  scattering  of  the  coherent  component,  does  not  lose  power  by  further  (multiple)  scattering  -  it 
scatters  into  itself  -  and  thus  dominates  at  large  distances  into  the  forest,  decreasing  at  a  much  reduced 
attenuation  rate.  In  the  transition  region  between  the  high  and  low  attenuation  regimes,  significant  beam 
broadening  and  pulse  broadening  occurs. 

In  this  study,  the  scalar  time-dependent  equation  of  radiative  transfer  is  used  to  develop  a  theory  for 
the  propagation  and  scattering  of  narrow  band1  pulsed  beam  waves  of  finite  cross-section  in  a  medium  that  is 
characterized  by  many  random  discrete  scatterers  (vegetation)  which  scatter  energy  strongly  in  the  forward 
scattering  direction.  Applications  include  the  scattering  of  millimeter-waves  in  vegetation  and  the  scattering  of 
optical  beams  in  the  atmosphere.  Strong  forward  scattering  occurs  at  millimeter  and  optical  frequencies  since 
all  scatter  objects  in  a  forest  or  in  the  atmosphere  are  large  compared  to  wavelength.  Again  of  interest  are  the 
range  and  directional  dependency  of  received  power,  pulse  broadening  and  distortion,  in  addition  to  the  effect 
of  a  finite  beamwidth  when  the  incident  field  is  not  a  plane  wave.  This  case  differs  basically  from  the  plane 
wave  case  in  that  scattering  out  of  the  beam  occurs  (while  in  the  plane  wave  case  ,any  multi  scattered  wave 
trains  will  always  remain  within  the  infinitely  wide  beam). 


i 


Narrow  band  refers  to  the  relatively  small  frequency  band  about  the  carrier  frequency,  which  permits  the  assumption  that 
media  parameters  are  independent  of  frequency. 
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2.  FORMULATION 


In  this  study,  vegetation  (a  forest)  is  modeled  as  a  statistically  homogeneous  half- space  of  randomly 
distributed  particles  that  scatter  and  absorb  electromagnetic  energy.  A  periodic  sequence  of  Gaussian  pulses  is 
taken  to  be  incident  from  free  space  onto  the  planar  boundary  of  the  forest.  The  incident  pulse  train  is 
assumed  to  be  either  (1)  a  cylindrically  collimated  or  (2)  a  spherically  diverging  beam  wave  (see  Figure  2.1). 

Chang  and  Ishimaru  [7]  used  scalar  transport  theory  to  study  the  scattering  of  a  monochromatic 
collimated  beam-wave  in  a  random  medium.  Their  approach,  however,  is  computational  intensive,  and  does 
not  provide  numerical  data  for  off-axis  beam  scattering1.  In  the  method  presented  here  -  which  also  involves 
using  the  scalar  transport  equation  in  cylindrical  coordinates  -  a  more  analytical  development  is  achieved, 
which  permits  numerical  data  to  be  obtained  for  off-axis  beam  scattering  in  the  forest. 

2.1  Incident  Gaussian  Beam  Wave  Pulse  Train 

For  the  cylindrically  collimated  case,  a  beam  wave  pulse  train  is  assumed  to  be  incident  from  the  air 
region  ( z  <  0 )  onto  the  planar  boundary  of  the  random  scattering  medium  (a  forest),  which  occupies  the  half¬ 
space  region  Z  >  0  (see  Figures  2.1  and  2.2).  At  the  origin,  the  magnitude  of  the  instantaneous  Poynting 
vector  of  the  incident  beam  wave  that  flows  in  the  +z  direction  is  assumed  to  be  given  by 

S(t)  =  2Savef  {t)  COS2  (*yc/) ;  f(t)  is  a  positive,  even  function  of  time  t  that  is  periodic  with  period  T, 
i.e.,  f(t  +  T)=  f  (t)  ,  and  normalized  such  that 


(2.1) 


Save  is  the  magnitude  of  the  time  average  incident  Poynting  vector  assuming  T  »  Tc  =  2 n  /  (Oc  ,  where 
COc  is  the  carrier  frequency. 

For  Gaussian  incident  pulses,  f(t)  takes  the  form 


(2.2) 


a0  = const . 


Limited  off-axis  data,  however,  is  available  in  [8]. 
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Since  the  incident  beam  wave  pulses  are  even,  this  periodic  function  of  time  can  be  represented  by  an  even 
Fourier  series  at  z  =  0 : 


/(<)=t+ 


(2.3) 


where 


2  71 


T 


2  772 

»  K  =—  /  f(t)cos(ixDt)dt 

1  -r  /  O 


(2.3a) 


Hence,  for  the  Gaussian  pulses  in  (2.2), 


fv  -  £VbV  -  s  e~{*vlao  f 


(2.4) 


a  o  has  to  be  chosen  properly  to  ensure  that  the  Gaussian  function  in  (2.2)  approaches  zero  as  t  —>  +  T/2 


allowing  the  limit  of  the  integration  in  (2.3a)  to  be  replaced  by  ±  oo  . 

The  specific  intensity  I inc  (power  per  unit  area  and  per  unit  solid  angle)  of  an  incident  collimated 

beam  wave  pulse  train,  that  travels  through  air  at  the  speed  of  light  "c"  in  the  positive  z  direction,  is  assumed 
to  have  a  p  -dependent  Gaussian  amplitude  with  beam  width  w  and  to  carry  power  density 


(2.5) 


with 


/((-z/c)  =  Re|X/„e'“'<‘-!/'j  (2.5.) 

and  fv  given  by  (2.  4)  for  the  Gaussian  pulses  (2.  2).  In  (2.5),  S(0)  is  the  Dirac  delta  function  and  6  is 

defined  as  the  scatter  angle  measured  positive  from  the  positive  direction  of  the  z-axis  (see  Figure  2.2). 

For  the  diverging  spherical  case,  the  specific  intensity  of  an  incident  spherical  beam  wave  pulse  train 
in  air  is  given  by 


(2.6) 
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with  time-periodic  pulse  train  f(t  —  R/c)  having  pulse  repetition  rate  (frequency)  MT  and  represented  by 
the  Fourier  series  expansion 


(2.6a) 


with 


R  =  [p2  +{z-z0)2]vl  ,  cos<90  =(z-z0)//?,  y/=(f)-(f)x  (2.6b) 

m)  is  the  directive  gain  of  a  transmitting  antenna  located  external  to  the  forest  in  free  space  on  the  z- 


(2.6b) 


axis  at  Z  =  z0  <  0 ;  Pave  is  the  average  radiated  power.  Spherical  coordinate  variables  ( /?,  dQ )  are  defined 

in  Figure  2.1b;  the  cylindrical  angles  ( (f),  (f>x  )  and  the  rotating  azimuthal  coordinate  angle  if/  =  (f)  —  (j>x  are 
defined  in  Figure  2.2. 

2.2  Phase  Function 

In  tranport  theory,  the  random  scatter  medium  is  characterized  by  an  absorption  cross-section  per  unit 
volume  <Ja  ,  the  scattering  cross-section  per  unit  volume  O  s  and  a  power  scatter  or  phase  function  p(s,s') . 

The  phase  function  depends  on  both  the  incident  power  unit  vector  direction  s'  or,  equivalently,  the  in-scatter 
angles  ( 0\ (/)' )  and  the  scatter  power  unit  vector  direction  s  or,  equivalently,  the  out-scatter  angles  (&,</>) 
for  each  scatter  event  (see  Figure  2.2). 

Since  scattering  surfaces  in  a  forest  have  essentially  random  orientations,  it  is  reasonable  to  assume 
that  a  forest  scatters  energy  symmetrically  about  the  direction  of  the  incident  radiation,  i.e.,  that  the  scattering 
which  occurs  at  each  point  in  a  forest  can  be  characterized  by  a  scatter  function  that  depends  on  the  angle  y 

between  s  and  s  ,  where  y  =  cos~ 1  (£'  •  i),  and,  therefore,  that  the  phase  function  can  be  written  as 


p(s,s')  =  p(s-s’)  =  p(y) 


(2.7) 


In  addition,  since  all  scatter  objects  in  a  forest  are  large  compared  to  the  wavelength  at  millimeter-wave  and 
optical  frequencies,  the  forest  scatters  energy  strongly  in  the  forward  direction  but  weakly  in  all  other 
directions.  For  that  reason,  the  scatter  function  is  assumed  to  be  characterized  by  a  strong  narrow  lobe 
superimposed  over  an  isotropic  background.  This  type  of  scatter  function  can  analytically  be  expressed  as  a 
Gaussian  function  added  to  a  homogeneous  term,  i.e.. 
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p{y)  =  aq{y)+^-a)  ,  q(y)  = 


—  ,Ays«x, 

AM 


(2.8) 


which  is  normalized  such  that 


jj  p(y)dQ.  =  4  n 


(2.8a) 


(K1  is  the  differential  solid  angle  about  the  scatter  angle  s  ,  A ys  denotes  the  width  of  the  forward  lobe  in 


the  scatter  pattern  and  a  is  the  ratio  of  the  forward  scattered  power  to  the  total  scattered  power.  The  scatter 
function  (2.8)  was  justified  in  [2]  by  reference  to  the  theoretical  and  experimental  comparison  of  results  in  [3, 
4],  and  by  the  experiments  conducted  by  Ulaby  et  al.  in  [9]. 

Therefore,  in  transport  theory,  according  to  the  above  discussion  and  taking  the  scatter  function  to  be 
specified  by  (2.8),  the  random  scatter  medium,  assumed  here  to  be  statistically  homogeneous,  is  characterized 

by  four  spatially  constant  parameters,  namely,  <Ja , (75,A/s  and  CC  .  These  four  parameters  are  understood  to 

be  "global"  parameters  in  that  they  remain  valid  at  all  points  in  the  random  medium  and  apply  to  an  average 
scatter  event  that  occurs  at  every  point  in  the  scatter  medium. 

2.3  The  Scalar  Time-Dependent  Transport  Equation  in  Cylindrical  Coordinates 

In  transport  theory,  the  specific  intensity  "  /  "  of  the  field  in  a  random  medium  is  governed  by  the  radiative 
transfer  equation  (transport  equation).  In  the  normalized  cylindrical  coordinate  system  ( p\z\l//}  for 
symmetric  scattering  about  the  direction  of  the  incident  radiation,  the  scalar  transport  equation  takes  the  form 


[10]: 


where 


V'  =  cos# —  +  sin#  cosy/ - sin#siny/ - 

dz  dp'  p'  dt// 

s-s'  =  cosy  =  cos  0  cos  #'  +  sin#  sin  #'cos(y/-y/') 


(2.9b) 


(2.9a) 


y/  =  (f) - (f)x,  dD!  =  sin#'d#'<iy/,' 


(2.9c) 
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The  normalized  space  and  time  variable  in  (2.9)  are  given  by  p'  =  <Jtp  ,z'  =  crtz ,  ty—Gtct  and 


W0  —  a sj  a t  ,  O',  =<Ja+<Js.  The  parameter  JTq  is  called  the  albedo  and  the  parameters  <Ja  ,  <JS  and  Gt 

are  the  absorption,  scatter  and  extinction  cross-sections  per  unit  volume,  respectively.  Implicit  in  writing  (2.9) 
is  the  assumption  that  all  parameters  are  independent  of  frequency. 

To  obtain  a  unique  solution  to  (2.9),  with  intensity  /  assumed  to  be  time  periodic,  requires  satisfaction  of 
the  following  two  boundary  conditions: 


/  =  Iinc  at  z'  =  0,  0<0<7i!2 
I  — >  0  as  z  — >  oo 


(2.10) 
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3.  RIGOROUS  SOLUTION  FOR  COLLIMATED  BEAM  WAVE 

As  is  customary  [1],  the  specific  intensity  is  separated  into  two  components,  namely,  the  reduced 
incident  intensity  I ri  and  the  diffuse  intensity  I d  by  letting 

I  =  ^ri+Id.  (3.1) 

Substituting  (3.1)  into  (2.9)  and  (2. 10)  yields  the  defining  equations  for  /n  and  I j ,  which  take  the  forms 


— Iri  +  s-V'Irt+Iri  =  0, 


(3.2a) 


— id  +  s  ■  via + /,  =  JJ  Pirv.dn  ■+  JJ  p(r)ijn  • .  (3.2b) 

with  boundary  conditions 


4* 


I  ri  ~  line  *  ^d  ~  ® 

/„.  ->  0 


7t 


at  z'  =  0,  O<0<- 

2  . 


(3.2c) 


as  z  — >oo 


where  s  •  V'  is  defined  in  (2.9a). 

To  solve  (3.2),  Fourier  series  representations  are  introduced  for  the  time  dependence  of  the 
intensities: 

Ij{pr,z\t\0,Y)  =  Re  j| ;4u(p',z',^^)e‘w' j  j  =  inc,ri,d  (3.5) 

where  Tf  =  <JtcT ,  o)9  and  normalized  angular  frequency  CO  —  (o/{c(7t) .  Note  that  in  the 

cylindrical  case,  both  Iinc  l)  and,  therefore,  /n  t;  depend  only  on  the  three  variables  p\z* ,6 ,  not  four 
variables  as  written  in  (3.5).  Observe  also  that  although  the  specific  intensity  (power  quantity)  is  always 
positive,  the  individual  Fourier  constituents  I j  v  may  be  negative  and  thus  cannot  physically  represent  power. 
Substituting  (3.5  )  into  (3.2)  yields  for  z’>  0 

(3.6a) 


,t><y  4i,t> + 5  •  v  i riv + iriv  =  o 


iuoi'h,v  +  *  •  +  L,o  =  P(yVd,u  sin  e  'dffdyf  ’+  ^  JJ  p(r)iri,u  sin  O' dG'dy/\ 3.6b) 


w. 


with  boundary  conditions 

4,=4«,  .  4*=°  <»  z’  =  °  >  O£0s| 

Iri  V  — >  0  ,  Id  y  ->  0  05  >  00 

where  L>  =  0,1,2,*  •• 

A  comparison  of  (2.5)  expressed  in  normalized  variables  with  (3.5)  gives 


(3.6c) 


j  _  _2  r  ^ ave  ~(p'/w')2  -iixo’z’ / cos 0 

1  inc,u  ~~  J u  ,2  e  o  •  n 

—  27rsin6 


7TW 


(3.7) 


where  S(0)  is  the  Dirac  delta  function  and  W  =  <Ttw  is  the  normalized  beamwidth;  note  that  in  (3.7), 
COS  0  =  1  because  0  =  0  due  to  the  Dirac  delta  function. 

The  first  boundary  condition  in  (3.6c)  together  with  (3.7)  and  the  definition  of  Iri  v  in  (3.5)  dictate 

that  lri  v(p\z\0)  oc  S(0)/[27rsin  0] .  It  then  follows  that 

sv'*ri*  =  cos0T-;L.o  +  sin  G  cos  y/  I ri - rsin<9siny/— / 

oz  op  p  dif/ 


a  d  7 

=  C°S0-/„.„ 


(3.8) 


and  (3.6a)  reduces  to 


(3.9) 


ivv’irtM +cos6>^7  =  o 

(3.9)  is  solved  by  direct  integration  and  use  of  boundary  conditions  in  (3.6c)  to  give  for  z'  >  0 

/  —rpf  -221  r~(p'/w')2 
ri,  v  ^ tJv  ,2  e  /,  •  a 

7tw  27ts\n6 

where  Pme  is  time-average  power,  p'  =  (Jtp  and z  =  a ,z ;  the  un-normalized  coordinate  variables  (p  ,6) 
are  depicted  in  Figure  2.2. 


(3.10) 


Let 


Ij,o(p'’z'>&>V)  =  IjJ,{p\z\0>V')e  va>z  ,  j  =  inc,ri,d 


(3.11) 


Substitution  of  the  reduced  incident  intensity  from  (3.10)  into  (3.6b,c)  and  using  (3.1 1)  gives  for  zf  >  0 


3  \  d  3 

[1  +  ivco'(\  -  cos  0)]Id  u  +  [sin  9  cos  if/ - sin  G  sin  if/ - 1-  cos  9  — ]/rf 

dp'  p’  dif/  dz' 

4  K  *'  4  K  7TW 


(3.12) 


with  boundary  conditions 


o 

II 

a 

at 

z'  =  0 

o 

as 

z'  — >00 

0<6  <  — 
2 


(3.12a) 


where  o  =  0,1,2," 

To  solve  (3.12),  Idv  is  represented  in  terms  of  the  Fourier- Hankel  transform 


00  00 

Id.o(p'’z',0’y/)  =  Yj  j  A^{k',z',6)\jm{k'p')cos(mif/)\k'dk'  (3.13a) 

m=0  kf=0 

with  inverse  transform 

2  n  oo 

A^{k\z\0)=  J  J  ldv(p\z\0i\f/)\jm(k'p')cos(mf/)\p'dp,dif/  (3.13b) 

^=0  p= 0 

The  representation  in  (3.13)  for  IdJU  is  an  expansion  in  terms  of  basis  functions  cos  m(k'p')  which 

are  complete  and  obey  well-known  orthogonality  properties.  The  0-dependent  coefficients  A” 
are  further  expanded  in  terms  of  Associate  Legendre  functions 

AZ(k\z\e)=fj£J(k\z')pr(™e)=fd(2i+\)£J(k\z')pr(™se)  o.m) 

l=m  l-m 

so  that  (3.13a)  becomes 

m=0  £'=()  l=m 

Since  scattering  is  assumed  to  be  symmetric  about  the  direction  of  the  incident  wave,  the  phase  function  is 
a  function  of  y  only  and  is  conveniently  represented  as  a  series  of  Legendre  polynomials  /}  [1] 

00 

p(r)  =  H(2l+x)sipt  (cos/)  (3.16) 

1=0 
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The  Legendre  polynomials  are  then  expressed  in  terms  of  Associated  Legendre  functions  via  the  expansion 

[1] 


n= I  (/  +  ft)  • 

= Ztt rT3T/r(//)/r(//')cos(/i(v' -¥')) 


(3.17) 


tieAl  +  nV- 


with  fJ.  =  COS  0  and  En  = 


f  2,  n  =  0 
1 1,  n  =1,2,... 


(3.17a) 


to  give 


1=0  n= 0 


en(i+ny: 


(3.18) 


For  the  scatter  function  (2.8), 
2a  r 


Si  A  2  I 

A  y  f 


°L  |  e^,Ar‘)2p,(cosy)smydr  +  (l-a)Sol ,  50l  =  j*’ 

ys  r= 0  l0’ 


1,  for  1  =  0 


(3.19) 


fori*  Oj 

Substituting  (3.15),  (3.18)  and  (3.19)  into  (3.12),  using  orthogonal  properties,  recursion  relations  and 
various  identities  of  the  functions^  (&'//) ,  COS  (my/)  and  Ff"  ( COS  #),  and  truncating  at  (m  =N  , 
l  =  N)  yields  the  following  inhomogeneous  system  of  linear  first  order  differential  equations: 


(l-nt) 


+(/  +  m  + 1) 


+(2/+i)[i-jp0gl+*w']<, 


+1[(/_OT_1)(/_OT) ^+u_,  -(/ + m + 1)(/ + m + 2) jCw ]  +  y  [K-u-i  ~ <-u+\ ] 


r-ju 


(3.20) 


=i^C/^*v/J),(2/+1)gA0^' 


for  m  =  0,  1,  2,  •••,7V  l  =  m,  m+\,  •••,7V  0<&'<oo 


with 


*«-.  =  < 


2,  m  =  1 
1,  m  =  2,3,. 
0,  m  =  0 


(3.20a) 


Introducing  the  normalization 


where 


2  1 

(/+m)! 

y  2/+1 1 

(/-*)! 

y=4^~m 


sm  = 


[2,  m  =  0 
1 1,  m  =  1,2,... 


(3.21a) 


allows  (3.15)  to  be  rewritten  as 


dju 


(p',s,o,v)=i }  <3-22> 


m=0  jt'-0  /=m 

and  simplifies  (3.20),  which  becomes 


a, 


— 7  Ki  l  -  iV0)bv  .  . 

m,/-l  m,l- 1 


+  a. 


+  »  f«i[aAtU-l  a5^m+l,/+l  ]  +  ^  a7^m-l,/+l  ]  —  Sl^mOe 


—b:M-ivco'b:J+i 


k' 


+  a3  [l  -  W0g,  +  iva)']b", 


(3.23) 


with 


£m  = 


dm  =< 


,  *  *  • ,  2V ;  l  =  m,  m  + 1,  • 

-,N; 

0<  k',z'  <00 

0  , 

O 

II 

5 

"I  =  °  <5 

m  =  l,2,3,...,JV 

F  , 
1  , 

ro  =  l 

m  =  2,3,...,N 

,  m  =  0 

,  m  =  1 

m  =  2,3,...,N -l,N 

II 

O 

Jl  ,  m=0 
[0  ,  m*  0 

(3. 23. a) 

a4  = 


«6  = 


+»X 

iTa) 

v 

(2/  +  1)2 

1  (/  -  m  • 

-lX/-w 

X2/-1) 

(2/  +  1)2 

l(/  +  m) 

(/  +  m  - 1 

X^TJ 

«X2/-l)  _  _  l(/-m  +  lX/  +  /H  +  lX2/  +  3)  _  _  /^T-T 

“75 - >a2  =  , - 7 - v> - ,a3  =  V2/  +  1  , 

1)2  V  (2/+1)2 


(2/  +  1)2 


,a5  =. 


a7  =  . 


(2/  +  1)2 


(2/  +  1)2 


(3.23.b) 
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3.1  Homogeneous  Solution 


The  solution  of  (3.23)  requires  determination  of  both  the  homogeneous  and  particular  solutions.  For 
the  former,  the  right  hand  side  of  (3.23)  is  set  equal  to  zero  and  the  homogeneous  solution  is  assumed  to  be  of 
the  form 

W>*')  =  <C  (3.24) 


Substitution  of  (3.24)  into  (3.23)  with  the  right  hand  side  set  to  zero  yields  the  homogeneous  system  of  linear 
equations 


for  m  =  0,  1,  2,  •••,7V’  l  =  my  m+ 1,  •••,7Vr  0<k'<oo 

(3.25) 

where  b ^  =  X-W^g/  +iuct)'  ,cv  =  iuco'  ;  (€m  ,Sm )  and an ,  n  =  1,2,..., 7,  are  given  in  (3.23.a,b), 
respectively.  Writing  (3.25)  in  matrix  form  gives  the  generalized  eigenvalue  equation 

[A0]G  =ct[C0]G  (3.26) 

where  matrices  [A0]  and  [Cq]  are  matrices.  The  eigenvalues  ( <J )  and  eigenvectors  ( G  )  are  determined 
by  using  the  QZ  method  algorithm  in  die  LAPACK  library  [11]. 

3.2  Particular  Solution 

The  particular  solution  to  (3.23)  is  obtained  by  assuming 

,3.27) 

which,  when  substituted  into  (3.23),  gives  the  system  of  linear  equations 
a\  +a2  +  aiPvl^m,l 

+yf.[“.Cu-l  +  f  <3-28) 

=  g,Sm0  for  m  =  0,1,..., N  and  /  =  m,m  +  l,...,N 
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where av  =14- ivco'  with  {^m^my8m0,bvl,  g{  and  an ,  n  =  1,2,..., 7  )  defined  in  (3.23a,b).  In  matrix 
form  (3.28)  is  written  as 

[b0  ]F  =  g  (3.29) 

where  [  B0  ]  is  a  matrix  while  F  and  g  are  column  vectors.  The  ZGESVX  routine  from  the  LAPACK 
library  is  used  to  find  the  particular  solution  [11]. 

The  general  solution  to  (3.23)  is  the  superposition  of  the  particular  solution  and  A^=(A^+l)  /4 

allowable  homogeneous  solutions  that  obey  the  condition  that  Re|l/a'|  >  0  ,  which  ensure  that  solutions 
decay  as  z  — >  oo  .  Thus,  the  general  solution  is  obtained  as 

+  £‘‘M(k,,z')  =  KAk') 

i=\  /=1 


e  *  -l- 


'L“iGZ{k')e-2',°‘  (3.30) 


where  an  i  =  1,...,  Nh  ,  are  arbitrary  constants  that  are  determined  from  the  boundary  condition  at  zf  =  0  . 
Thus,  using  (3.30)  in  ( 3.22)  and  incorporating  the  truncations,  the  diffuse  intensity  is  expressed  as 


m=0*'=o  /=m 

/r(cos^)cos(myr) 

u;  Vm 

3.3  Boundary  Condition  at  Entrance  to  Forest  (z'  =  0)  to  Find  the  Coefficients  a. 

From  (3.12a),  Id  v  must  satisfy  the  boundary  condition 


(3.31) 


Idu(p\z'  =  0,0, ^)  =  0  for  z  —  0,  O<0<|,  0  <  p'  <  oo,  0  <  y/  <  In  (3.32) 


Introducing  (3.31)  into  (3.32)  gives 


Idu{p',z'  =  0,0,y,)  =  J  Idu (k',z'  =  0,6,y/)Jm  (k'p')k'dk'  =  0 


k'=0 


(3.33) 


where 


w.vVL  X  (2/+r )ff(k\w)b:,(k\z') 


m= 0  l-m 


/f  (cos#)  cos(m^) 

uf  Vm 


(3.33a) 


with 
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(3.33b) 


C(*^>Cr+I>,<W 


1=1 


Muliplying  (3.33)  by  J m{k"  p') ,  integrating  over  p'  from  zero  to  infinity,  interchanging  the  order  of 
integration  and  using  the  completeness  relation 


j  J. (k'p')Jm  (*' V) p'dp'  =  6(k'-k’)lk' 

P-  0 


(334) 


yields 


au 

Id.»{k'’z'  =  °’0’'/')=  J  IdAp'^’  =  p')p'dp'  =  0  (3.35) 


p'=0 


Thus  with  (3.33a),  the  following  boundary  condition  is  obtained: 


m=0  /=m 


(336) 


;r 


forz'  =  0,  O<0<— ,  0</rf<oo,  0<i//<2  n 


The  boundary  condition  (3.36)  is  satisfied  by  using  the  normalized  spherical  harmonic  functions  as 
testing  functions  in  the  weighted  residual  method  [12].  Hence, 


2  n  n/2 

J  }  J**  (*'.*' =  ^0.^) 

^=0  0=0 


~V(oos0) 

cos  ( my/ ) 

U? 

J/™ 

sin  QdGdif/  =  0  (337) 


for  values  of  (/  —  m)  which  are  odd.  Substituting  (3.33a,b)  into  (3.37)  yields 


N  V 


e"(t'.«,)ZZ(2',+1) 


m'=0  r=m' 


Kf+'LvpJ, 


/=! 


2k 

\ 

y/=0  _ 


cos  (mV) 


cos 


T 

~/f'(cos0)‘ 

/f  (cos#) 

sin  0d6 

(9=0 

uf 

u; 

(3.38) 


•S'  = 


Thus, 
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!«■(*') 


/=! 


S(2/'+i  )<%,(*")  fr 

J'=m 


-’Z(2l'  +  l)F^(k')/~ 


V—m 


•W) 


or,  equivalently, 


where 


1=1 


«:;(*')= z  (a- +i)g" 

V-m 

1 

II 

1? 

with 

-  "f 

7,7=1 

(9=0 

Fp  (cos  9) 

up 

pp  (cos#) 

up 

f  A 

yv 

SI 

V-m 


sin  GdO 


J”j.  for  m  >  0  such  that  l-m  =  odd, l - 1'  =  odd 
y  form>0  such  that  l-m=odd,l-l'  =  0 
0  otherwise 


(3.39) 


jm  =  J_  ^  j  y r+/)/2+m+ 1/2  V 


(2/  +  l)(2/'  +  l)^(/  +  /n)2  -m2  j 


(r-/)(r+/+i) 

R(r-m)R(r  +  m)R(l  +  m  +  l)R(l-m  +  \) 


x=even 


(3.39a) 


Assume  that  a  highly  directive,  lossless  antenna  of  narrow  beamwidth  and  narrow  bandwidth  is 
located  inside  the  forest  [2].  This  receiving  antenna  is  characterized  by  an  effective  aperture  A(yR),  where 

Y R  is  the  angle  included  between  the  direction  of  observation (0R,y/  R)  and  the  pointing  direction  of  the 
antenna  axis,  i.e.,  the  main  beam  direction  (Qm  W m  )  ;  see  Figures  2.2  and  2.3.  Hence, 


3.4  Power  Received  by  a  Highly  Directive  Antenna  in  the  Forest 

cos yR  =  cos0Rcos0M  +sm0Rsin0M  cos (y/R  (3.40) 

In  transport  theory,  powers  add.  Thus,  the  instantaneous  power  received  by  the  antenna  is  the  sum  of 
the  intensity  contributions  coming  from  all  directions  multiplied  by  the  effective  aperture  of  the  antenna,  i.e., 
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p,(z',p',t',eu,¥u)=^\f1p;(,z',p',eu,Vu)e‘"v-^ 

1^=0 


(3.41) 


where 


(z  *P  Wm)  ~  \\  Ae{y r)I , p  ^0Ryy/R)s\n dRdORdy/  R  (3.41a) 

4  x 


and 


K  ”  Idtv  +  Ifi%v 


(3.41b) 


Note  that  6  —6 R  and  y/  =y/R  . 

For  millimeter  waves,  the  carrier  frequency  is  very  large  and,  therefore,  the  bandwidth  of  the 
received  signal  is  narrow.  For  such  a  small  bandwidth,  the  effective  aperture  and  gain  of  the  receiving  antenna 
can  be  taken  to  be  independent  of  frequency  and  to  be  related  by  the  general  expression 


Ae(rK)  =  ^D(yR) 
4  71 


(3.42) 


where  XQ  is  the  free  space  wavelength  and  D(y R)  is  the  directive  gain  of  the  antenna  at  the  carrier 
frequency. 

Assume  that  the  directive  gain  is  Gaussian  with  a  narrow  beamwidth  and  no  sidelobes,  i.e.. 


f  o  V 


D{yR)  = 


\^y  m  j 


,  A  yM  «  K , 


(3-43) 


which  is  normalized  such  that 


j\D(rR)sindRd9Rdi/R  =4n 

4  K 


(3.43a) 


Using  (3.41b),  (3.42)  and  (3.43),  the  received  instantaneous  power  PR  in  (3.41)  is  obtained  as  the 
sum  of  diffuse  power  Pf  and  reduced  incident  power  P*  .  The  received  diffuse  power  is  obtained  as 

P;(z\p\t\eu,¥u)  =  Ke{fiPf„(z\p\eu,¥„)el'*'-^\  (3.44) 


0 


where 
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(3.44a) 


pd*v(z'’  P\  0M,yM)  =  \\  MYr  )Id,u(z\  p\  Qr,Vr)  sin  ORd0Rdy/R 

4  n 

=  ^  JJ  °(Yr  )IdAz',  p\  Or  Jr)  sin  0RdORdy/R 
£ 


= -f- !d.u  (z'>  p\  0M ,  VM )  jj  D(rR)  sin  0Rd0Rdy/R 

4  n 

=  \^d,v(z  *P  • 

Similarly,  the  received  reduced  incident  power  is  obtained  as 

P‘(z\p\t\eu,Vu)  =  Kt\fiP^v{z\p\eu,Vu)e‘"v-1’A  (3.45) 


L^=o 


where 


>P  —  \\ ^e(y r)I * p  iQR,$R)sinORd9Rd\i/ R 


4  n 

erf 

Pa. 

An 

7lW2 

Jo 

ii 

<7, 2 

Pa. 

An 

nW2 

Jo 

In  n 


2;rsin0o 


The  instantaneous  received  power  is  normalized  to  the  received  time-averaged  power  at  z'  =  0 ,  p'  =  0  , 
0M  =  0  and  —  0 ,  which  is  given  by 

(pr (0, 0, t 0, 0)) = T  J  pR (0,0,t’,0,0)dt'  =  -^-D(0)o? 

1  _f/2  n  7rv 


_  ,  (3.46) 

4n  kw 


since  Id  v  is  zero  at  z'  =  0  and  using  (2.1)  and  (2.3).  Thus,  the  normalized  total  instantaneous  power  is  the 
sum  of  the  reduced  incident  and  the  diffuse  normalized  received  powers,  namely. 


P\z'  p' t' 0  w  )=  ^-SLlRA-  '0mW_m1  =  p>  +p> 

' ,p’ • u'¥u)  (p"  (0,0,1' AO))  *  -■ 


(3.47) 


Using  (3.44-3.47),  the  total  normalized  instantaneous  received  power  becomes 


where 


pI(z',p',«„,Vh  ) = c + f»e'' + ) 
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(3.48a) 


Observe  from  (3.48)  that 


K 


u= 0 


e-(p'/w’)2 c~z' 

D(  0) 


v—0 


(3.49) 

£>(0) 


This  result  was  expected  since  the  reduced  incident  power  P\  is  the  power  of  the  incident  beam  wave  pulse 
train  which,  as  it  travels  through  the  forest  along  a  straight  path,  decays  exponentially  due  to  absorption  and 
scattering,  but  maintains  its  narrow  beamwidth.  The  6  -dependence  of  thus  reproduces  the  radiation 

pattern  of  the  receiving  antenna.  On  the  other  hand,  the  diffuse  intensity,  generated  by  scattering  of  the 
reduced  incident  intensity  and  by  self-regeneration  due  to  multi  scattering,  is  characterized  by  a  broad 
beamwidth  which  is  larger  than  that  of  the  receiving  antenna.  Hence,  the  receiving  antenna  acts  to  probe  the 

angular  distribution  of  the  diffuse  intensity  Id  u  as  seen  from  (3.48).  Note  that  at  z'  =  0 ,  p'  =  0  and  for 

9 M  =  0  ,Pri  =  f  (/')  .  In  numerical  results,  the  summation  in  (3.48)  is  truncated  at  a  value  O  =  Umax  that 
ensures  convergence. 
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4.  RIGOROUS  SOLUTION  FOR  SPHERICAL  BEAM  WAVE 


The  important  practical  case  of  a  spherical  wave  (radiated  by  an  antenna  placed  in  the  air  region 
external  to  the  forest  half-space)  is  assumed  to  be  incident  upon  the  forest.  Such  an  incident  beam  is  given  by 
(2.6)  and  illustrated  in  Figure  2.1b.  In  (2.6),  the  antenna  power  pattern  or  directive  gain  D  is  assumed  to  be 
rotationally  symmetric  so  that  it  depends  only  on  the  angle  0Q . 

Following  the  development  of  the  rigorous  solution  for  the  collimated  beam,  the  specific  intensity 
for  the  spherical  beam  is  assumed  to  be  the  sum  (3.1)  of  the  reduced  incident  intensity,  defined  by  (3.2a),  and 
the  diffuse  intensity,  which  is  governed  by  (3.2b),  that  satisfy  the  boundary  conditions  (3.2c). 

Introducing  the  Fourier  series  representation  (3.5)  into  (3.2)  yields  the  defining  equations  (3.6)  for 

the  transformed  intensities  Ij  o(p\z\0,y/\  j  =  inc,ri,d.  Using  the  expression  for  the  incident 
spherical  beam  (2.6)  with  the  Fourier  series  representation  (3.5)  gives 

i  0)<%o  —iuccf (z'—Zq)/ cosO  .v 

AkR,2  sin0  (4.i) 


where  R'  —  <J(R ,  z'0  =  <Jtz0  with  ( /?,  #0 )  defined  in  (2.6b)  and  in  Figure  2.1b. 

Use  of  (4.1)  and  boundary  conditions  in  (3.6c)  give  the  reduced  incident  intensity  for  the  spherical 
beam  for  zf  >  0  as 

7  -2 ,  -( 1+iixw')  r'/cos  Oq  iuw'z'n  /  cos  6q  S(9-90)S(y,) 

InJv~utJu  A  n,2  e  e  n 

An R  sin  0O 

where  Pm  is  time-average  radiated  power  and  zQ  <  0  .  Equation  (4.2)  can  be  shown  to  satisfy  (3.6a). 

Substitution  of  the  reduced  incident  intensity  from  (4.2)  into  (3.6b,c)  and  using  (3.1 1)  gives  for 

z  >0 


^1  d  ^ 

[1  +  iuco'(l  -  cos  9)]Id  +  [sin  9  cosy / - sin  9 sin  yr - +  cos  9  — ]/, 

dp'  p'  dy/  dz'  ‘ 


An 


An ' 


AnR' 


(4.3) 


with  boundary  conditions 


20 


71 


Id,u=0  at  Z'=0  t 

Idv  — >  0  as  z'  —>  oo 


(4.3a) 


where  =  0,1,2,-  •• 

To  solve  (4.3),  Id  v  is  represented  as  in  (3.22)  in  terms  of  the  Fourier-Hankel  transform  and 
Associate  Legendre  functions  as  follows: 

)  i(2<  +  I)<r (*•■  -X, (f .*)■ K (JTa) (4.4) 

m=0  k'=0  l=m  l  ' 


where  U™  and  Vm  are  defined  in  (3.21a)  and 


W 

Q°  =  — p  f 
^  Anvl  aveJv 


(4.4a) 


Substituting  (4.4)  into  (4.3)  using  the  phase  function  (3.18)  with  (3.19)  gives  the  inhomogeneous 
linear  system  of  equations  that  needs  to  be  solved 


a, 


—buml\-iuG>'buml  i 

t  mj-l  m,l— 1 


+  a. 


^7  Km  -iv(0'K.M 


+  a3  [l  -  fV0g,  +  iu(o']bvml 


(4.5) 


k  -  r  -|  kf  r  -l 

+  ~Z£m  L**  A+1./-1  —  a5^m+\,l+\  J  +~Z^m  \_afPm-\J-\  ~  a7^m-l,/+l  J  =  7  ~7Z)  .  ix-il/2  Si  Pi  >  Z  ) 

2  2  K,(2/  +  1)J 


m  =  0,  1,  2,  •••, N ;  /  =  w,  m  + 1,  0<k',z'<oo 


with 


£•„  = 
m 


5  = 

m 


0  ,  m=  0 

>/2  ,  /n  =  l  =• 

1  ,  /w  =  2,3,...,  N 


\  yfl  ,  /n  =  0 

[l  ,  m=\,2,...,N 

2  ,  m  =  0 


ll  ,  w  =  l,2,...,A( 


= 


0 

1 


,  m  =  0 
,  m  =  1 

,  m  =  2,3,...,JV-l,Ar 


^m0  —  ‘ 


0 


m-0 
m*  0 


(4.5a) 
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and 


_  l(/-mX/  +  mX2/-l)  _  /(/-«  +  lX/  +  w  +  lX2/  +  3)  _  p— 

J — ^ —  ’“2 '  i — - ,aj  -  ^ + 1  • 


(2/ 


+  ] 


_  to-m-\ll-mpL-\)  _  I(/  +  /«  +  lX/  +  m  +  2X2/  +  3) 

a4=i~~^if  • “5i  (^1?  ■ 


«6  = 


/(/ -t- mX/  +  m - lX2/ - 1)  _  |(/-m  +  2X/-w  +  lX2/  +  3) 


(2/  +  1)2 


al  = 


(2/+O2 


(4.5b) 


/f  (COS0O) 


P,-(*',z')=  j  ''  V,7;"°v.(t>vv 


(4.5b) 


p’= 0 


Note  that  the  linear  system  (4.5)  is  identical  to  (3.23)  for  the  collimated  beam  except  for  the  forcing 
term  that  is  proportional  to  (4.5b).  The  homogeneous  solution  of  (4.5)  thus  remains  the  same  (and  is  given  in 
(3.30)) ,  but  the  particular  solution  differs. 


4. 1  Particular  Solution 


Although  the  spherical  beam  appears  very  similar  to  the  collimated  beam  in  the  formulation  of  the 
solution  ,  the  computational  complexity  involved  is  significaant.  Two  approachs  to  finding  the  particular 
solution  are  presented  below. 

4.1.1  The  Fourier-Cosine  Approach 

Let 

PmjW,z')=e-'GmJ(k\z')  (4.6) 

with 


J-n  K 


e-z'/cos00+z'  /y»( coS0o) 


Jm(k'p')p'dp' 


(4.6a) 


rrm 

p'= 0  MX  Ul 

In  (4.6),  the  exponential  factor  e~z  was  removed  from  the  integral  over  p'  in  the  forcing  term  p™  (k\ z)  in 
the  system  of  equations  (  4.5).  The  remaining  integral  (4.6a)  still  approaches  zero  as  z  — >  00 ,  i.e., 
Gm  z)  — >  0  for  z  — >  00  as  required  for  the  Fourier-Cosine  transform  to  exist. 


The  Fourier-Cosine  transform  of  Gm  t{k\ z)  takes  the  form 
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00  1  00 

G„,(*',z')=  J  G„,(k\u)cos(uz')du  =  -  \  GmJ(k’,u)el“'du  (4.7) 


M=0 


where 


GmJ(k',u)  =  Gml(k',-u)  =  —  f  G  j(k',z')cos(uz')dz' 

*/-o 

=  mL  L  vr~ 


Jm(k'p')cos(uz')p'dp'dz'  (4.7a) 


z'-O  p'=0 


Thus,  the  forcing  term  in  (4.5)  for  the  coefficients  b”  {  becomes 


p„J(k\z')  =  e-GmM’,z’)‘  \  } 


(4.8) 


In  this  form,  the  z-dependence  of  Pmj(k\z')  is  exponential,  as  desired,  and  the  calculation  of  the  particular 
solution  follows  standard  procedure,  i.e.,  one  assumes  that 


f>:;f (*’.*')-  j  b:Hk\u)e"-‘"'-du 


(4.9) 


And,  therefore. 


e  b:.p(k\z')  * 


mj 

dz' 


=  -  j  [1 


(4.10) 


Substitution  of  (4.9)  and  (4.10)  into  (4.5)  yields  the  following  linear  system  of  equations 

«■  [-P  - iuKj-i  ~ iva'K!- 1 ]  +  «2  [~P  - iuK.M  -ivco'KL \ ]  +  «3  P  “  wo8i  +  **»']*£f 


(4.11) 


+ 7  ]  +  4-  [>44&-,  -  ]  =  ,.-7 7,  g,G..,(k» 

2  L  2  L  J  2[£m(2/  +  l)J 

m  =  0,  1,  2,  /  =  m,  /w  +  1,  0<  k\z'  <00 


where  all  parameters  have  been  defined  in  (4.5). 

Once  the  linear  system  (4.11)  is  solved,  the  particular  solution  is  found  via  (4.9).  Because  of  the 
common  exponential  factor  e  z  in  (4.9),  the  particular  solution  decreases  with  z  as  required.  Partial  fields 
that  increase  with  z'(  which  in  principle  are  allowed  by  the  transport  equation  but  are  unphysical)  are 
automatically  excluded  in  this  formulation.  Since  the  exponential  factor  was  removed  from  the  modified 
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forcing  term  Gm  t{k\z)  ,  this  term  decreases  slowly  with  z  and  extends  over  a  relatively  wide  z*  -range. 

Hence,  its  Fourier  transform  Gml(k\u)  is  relatively  narrow-band  in  W,  which  reduces  the  numerical 

burden  associated  with  the  computation  of  the  overall  particular  solution  since  only  narrow  band  of  u  values 
are  needed.  Once  the  particular  solution  is  determined,  the  final  specific  intensity  and  the  received  power  are 
found  by  proceeding  in  the  manner  described  in  Chapter  3. 

The  method  presented  above  for  the  spherical  beam  is  computationally  demanding.  This  is  because 

calculation  of  Gm  l(k\u)  requires  evaluation  of  a  double  integral  over  p'  and  z  via  (4.7a)  and  the 

determination  of  the  overall  particular  solution  then  requires  an  additional  integration  over  li  (though  over  a 
limited  u  -range).  In  addition,  a  summation  over  V  is  needed  to  obtain  the  time-dependent  specific  intensity. 
In  this  report,  only  the  time-  independent  case  is  numerically  treated  for  the  spherical  beam. 

4.1.2  Finite  Difference  Solution 

An  alternative  approach  to  finding  the  particular  solution  for  the  spherical  beam  case  is  to  use  a  finite 
difference  scheme  to  solve  (4.5).  In  this  approach,  the  z  -coordinate  is  discretized  by  setting 
z'  —  jAz\j  =  0,1,2,... ,  and  the  differential  quotients  in  (4.5)  are  replaced  by  the  difference  quotients: 

a?  Z?  (412) 

To  solve  the  resultant  system  of  equations,  the  finite  difference  method  requires  values  for  all  the 
coefficients  in  a  boundary  plane  Z  —  Z .  This  is  provided  by  approximating  the  diffuse  intensity  in  the 
boundary  plane  using  first  order  multiple  scattering  theory.  This  scattering  theory  yields  good  results  for  the 
diffuse  intensity  at  small  penetration  depths  into  the  forest.  Hence,  the  boundary  plane  can  be  chosen  to  be 
located  at  z  =0  or  1,  for  example,  whichever  is  more  accurate. 

With  the  diffuse  intensity  Id  «  I™(p\z'  =  z\ 0,y/)  known  from  first  order  theory  (the  time- 
independent  case  is  only  considered,  which  means  that  V  =  0  in  (3.5)),  the  coefficients  b%p  (k\ £')  are 


found  by  using  the  representation 


where  Q°  (k' ,w'}b“j  (k' ,z')  =  Q°  (k\w')b„j  (k' \z)  for  o  =  0  and  then  using  the  orthogonal 
relations 


}  PT  (cos 0)  P?  (cos 0)  2f  cos ( my/ )  cos  ( mV ) 

|  - *-JLA - -sm0d6  =  S,r  ,  f  - - - — r—^-dy/  =  8  ■  , 

J  j  jm  j  jm  «»«  J  jrm  rym  •  m.m  ’ 

0=0  u/  ur 


^=0 


,_„r  S(k'-k') 


J  JJk'p')JJk'p')p’dp'  = 


p'= o 


(4.14) 


to  obtain 


JtJ  J  It°  (p',i',9,v) P>  a>Sj,T>'^sinffded,//iJm  (k'p’)p'dp'  sin  $dff  dyr 

p'= 0  ip=0  9=0  U  t  " 

(4.15) 

Again,  once  the  particular  solution  is  obtained,  the  procedure  in  Chapter  3  is  used  to  find  the  specific 
intensity  and  the  received  power. 
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5.  APPROXIMATE  SOLUTIONS  OF  THE  TIME-INDEPENDENT  TRANSPORT  EQUATION 

5.1  Introduction 

Although  very  accurate,  the  rigorous  solutions  presented  are  computationally  intensive  and  take 
considerable  computer  time  to  generate  meaningful  data.  They  are  thus  limited  in  their  utilization  for  real  time 


applications  as  are  often  needed  in  military  situations.  To  provide  such  real  time  data,  approximate  solutions 


are  invaluable  even  though  they  are  not  as  accurate  as  the  rigorous  solutions.  Three  approximate  solutions  are 
developed  here  for  both  the  collimated  and  spherical  beams.  The  first  one  is  the  first  order  multiple  scattering 
solution  which  is  valid  at  short  penetration  depths  into  the  forest  half-space.  The  second  approximate  solution 
is  an  asymptotic  solution  that  is  valid  at  large  penetration  depths.  The  third  approximate  solution,  which  is 
referred  to  as  the  composite  approximate  solution,  is  a  combination  of  the  first  two.  These  solutions  are 
simpler  than  the  rigorous  ones  and  are  capable  of  providing  real  time  results  for  the  problem  of  propagation 
and  scattering  of  millimeter- waves  in  a  forest  region. 

5.2  First  order  Multiple  Scattering  Solution 

The  time-independent  transport  equations  in  cylindrical  coordinates  for  the  reduced  incident  and 
diffuse  intensities  along  with  boundary  conditions  are  specified  in  (3.2).  For  first  order  scattering,  the  phase 
(scatter)  function  in  the  second  integral  on  the  right  side  of  (3.2b)  remains  the  same,  however,  the  phase 
function  in  the  first  integral  is  approximated  by  a  8  -function  in  the  forward  direction.  Physically,  the  second 
integral  depicts  coherent  scattering  into  incoherent  scattering,  which  is  first  order  scattering,  while  the  first 
integral  describes  incoherent  scattering  into  itself  which  is  higher  order  scattering.  The  justification  for  using 
the  delta  function  to  characterize  the  higher  order  scattering  in  the  transport  equation  is  based  on  the 
assumption  that  particles  (leaves,  twigs,  etc.)  in  the  forest  scatter  strongly  into  the  forward  direction.  Hence, 
the  phase  function  in  the  first  integral  on  the  right  side  of  (3.2b)  is  assumed  to  be  given  by 


(5.1) 


where  a  is  a  parameter  that  is  chosen  to  optimize  this  approximation. 


Substitution  of  (5.1.1)  into  (3.2b)  yields 
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(5.2) 


s-vid+gli1{p\^,z\e,^='^l\p(r)irl(p\^,z',eJ)dCi 


where 


<7,  =  l-aW0  (5.2a) 

cos/ =  cos#cos0+sin#sin0cos(^-^  )  (5.2b) 

and  d  /  d  \j/  =  —d  /  d(j>x  is  used  in  the  definition  of  the  operator  s  •  V'  in  (2.9a)  since  {//  =  (f)  —  (f>x . 

The  solution  Id°  to  (5.2)  that  satisfies  boundary  conditions  (3.2),  called  the  first  order  multiple 
scattering  solution,  can  be  shown  to  be  given  by 

i y  z'  {*-*') 

C(p',t„z',e,t)=  — >—  f  e  '  “*  [ff  p(r)l„(p',l,l'^)sin0d0dm' 
4xax0  " 


for  0  <  0  <  n  /  2  (5.3a) 

"0  J  ,0 j)sin0d0dj]dz' 


W  f  -a 

i;u(P',<f>x,z',e, <(>)=-  ^ 


An  cos  0 


4  n 


for  n/2<0<n  (5.3b) 


where  6,  <f>,  z  are  integration  variables,  y  is  defined  by  (5.2b)  and 


p'  =  [p'2  +  (z'-z')2  tan2  ^-2/?'(z'-z')tan^cos(^-^t)]12 


p'  sin  -  (z'  -  z')  tan  0  sin  <f> 
p'  cos  <j>x-(z'-z’)  tan#  cos  ^ 


(5.3c) 


(5.3d) 


Note  that  I  d°  is  continuous  at 0  =  n / 2 . 

5.2.1  Collimated  Beam 

For  the  collimated  beam  having  a  Gaussian  beam  profile,  the  reduced  incident  intensity  (3.10) 
with  V  =  0  is  given  for  z'  >  0  by 


iri=ir,0U>',~z',d)  =  o2 


Pave 

nw'2 


e^H)  e~r 


2;rsin# 


(5.4) 


Substitution  of  (5.4)  into  (5.3)  gives 
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I™(p',z',e,¥)=  <j2  Wf™p{e)  f  for  0  <  6  <  7Z /  2  (5.5a) 

An w  cos#  ~?=0 

IFO(p',z',0,y/)  =  -<j2  Wf™p(0)  f  e^-^'e^'^dz'  ioxnl2<6<n  (5.5b) 

An  w  cos  0  -iz, 

where  if/  =  (f)  —  (f>x  and  p  is  given  by  (5.3c). 

The  solution  (5.5)  is  a  first  order  multiple  scattering  solution,  but  with  an  improvement  in  that  higher 
order  scattering  is  included  because  of  the  approximation  (5.1)  which  is  justified  because  the  phase  function 
has  a  strong  forward  lobe. 

5.2.2  Spherical  Beam 

For  the  spherical  beam  having  a  rotationally  symmetric  pattern  function,  the  reduced  incident 
intensity  (4.2)  with  u  =  0  is  given  by 

p„,d0«)  s(e-ea)5(i-l) 


/„=4  = 


4n  R' 


sin  Gn 


(5.6) 


where 


cos 4  =(z'-z'0)/ R' ,  sin  4  =p' / R' ,  R'  =  [p'2  +(z  -z')2]1 

Upon  substitution  of  (5.6)  into  (5.3)  it  can  be  shown  that 


(5.6a) 


f Fd°(p',z\6,V)=  f  forOSOS^  (5.7a) 

(An)  cos  0  ~,_Q  R  2 


IFO (p\ z'\G,y/)=- erf  Wfave  f 
d  y  ‘  (4nf  cosO ;{,  R'2 


>  e  ''  cos 8  1  [m0°]dz'  for  y  <  6  <  n  (5.7b) 


where  If/  —  (f)  —  (px  and 


cos/ =  cos#cos#o+sin0sin#o  cosy/,  \j/  =(!>-<f>x  (5.7c) 

Note  that  in  the  integrands  in  (5.7a,b),  the  expressions  COS  y  ,  COS  0O ,  sin  60  and  R'  are  functions 
of  all  the  variables,  i.e.  they  are  functions  of  p\  z\ 0,lf/  and  of  the  integration  variable  z  . 
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5.3  Asymptotic  Solution 

5.3.1  Case  1 

In  the  forest  at  distant  observation  points  on  or  near  the  Z  -axis,  the  reduced  incident  intensity  is 

negligible  while  the  diffuse  intensity  decreases  exponentially  withz*,  is  independent  of  p  and  <f>xy  and 

exhibits  a  scatter  pattern  that  depends  only  on  0 ,  but  not  on  (f)  (observations  that  were  made  on  examination 

of  data  of  the  rigorous  solution  for  the  collimated  beam).  In  the  development  of  an  asymptotic  theory  for 

z  » 1,  p* « 1  in  the  range  0 <0<7T,  the  time-independent  transport  equation 

in  cylindrical 

coordinates  (3.12)  for  diffuse  intensity  1  d  =  Id  0,  V  =  0,  was  shown  to  reduce  to 

cos^|7/rf,0+/,>0=fjjp(r)/rf,sin^^' 

L  4  it 

(5.8) 

with  boundary  condition 

I d  q  — >  0  as  z  — >  qo 

(5.8a) 

where  COS  y  is  given  in  (2.9b). 

For  z'  »  1  and  p' «  1  ,  assume  that 

Id=Id^IZ^^zSJcosG) 

(5.9) 

with  scatter  pattern  S^COS#)  taken  to  be  a  function  of  cos  0  that  is  normalized  to  unity 

for  0  =  0  and 

where  is  a  positive,  real  constant.  Substitution  of  (5.9)  into  (5.8)  gives 

W  rr 

(1  -  cos  0)S n (cos 0)  =  —  [|  p{y)S n (cos  0f) sin  ffdffdy/' 

4  n  y 

4k 

(5.10) 

To  solve  (5.10),  it  is  convenient  to  represent  both  S^COsG)  and  p(y)  in  terms  of  Legendre 

polynomials.  Hence,  assume 

SJcosfl)  =  X(2/+iy>,/>(cos0) 

1=0 

(5.11) 

and,  as  in  (3.16),  assume 

p{y)  =  Y,{2l  +  \)g,Pl{cosy) 

1=0 

(5.12) 

where,  as  before  in  (3.19)  for  the  scatter  function  (2.8), 
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g,  =  |  e'(r/Ar' ,2  P, (cos  y)  sin  ydy  +  (  1  -  a)£0/  ,  £0/  =  ^  (5. 1 2a) 

*/,  r{ o  [0,  /or/*0j 

Using  the  expansion  (3.17)  to  convert  (5.12)  into  the  form  (3.18)  reduces  the  double  integral  in 

(5.10)  to  a  single  integral  over  O'  since 

1  2?  00 

—  f  Pir)dw'  =p(cos^,cos^)=E(2/+l)^(cos<9)^(cos6>')  (5.13) 

ln  /=0  /=o 

Hence,  (5.10)  becomes 

If  r 

(1  -  cos  9)SV  (cos  0)  =  —  [  p(cQs6,cosQ')Sx(cosO')s,\nO'dO'  (5.14) 

^  ^=o 

which  is  an  integral  equation  for  S^COS  6 )  . 

Substituting  (5.11)  and  the  summation  in  (5.13)  into  (5.14)  and  using  orthogonality  and  recursion 
relations  yield  the  linear  system  of  equations  for  the  expansion  coefficients  bl : 

lb,_x  — z-  (2/  + 1)(1  -  W0g! )  +  (/ + 1)6/+1  =0,  /  =  0,1,2,...  (5.15) 

^00 

The  eigensolution  bt  to  the  linear  system  (5.15)  associated  with  the  smallest  real  eigenvalue  min  is  used 

to  construct  the  solution  S^COS  6)  in  (5.9).  All  the  other  larger  real  eigenvalues  decay  faster  and  are  ignored 

in  the  far  zone.  Because  the  diffuse  intensity/^  physically  represents  power,  it  is  a  positive  quantity  and, 

therefore,  the  solution  S^COS#)  used  in  (5.9)  is  also  positive  for  all  0  <  6  <  71 .  As  is  shown  in  the 

numerical  results,  Sm  (COS  0)  is  characterized  by  a  strong  maximum  in  the  forward  direction  0  =  0  and  is 

small  in  the  backward  region  n  /  2<G  <  7Z  . 

5.3.2  Case  2 

A  second  asymptotic  solution  to  the  time-independent  transport  equation  is  developed  below  for 
z  »  1,  p ' »  1  with  the  requirement  that  tan  0O  =  p'/(z'-z'Q  )  remain  sufficiently  small  in  comparison  to 

the  angle  width  of  (cos#);  z'Q  and  0O  are  defined  in  Figure  5.1.  This  latter  constraint  confines  the 
observation  point  in  a  conical  region  about  the  beam  axis.  Hence,  the  diffuse  intensity  is  assumed  to  be 
symmetric  about  the  direction  0O  rather  than  0  =  0  and  to  take  the  form 
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Id  =  Id.o  *  IZ  s  )  (5.16) 

withS^COS#,)  normalized  to  unity  for  9X  =0and  where  ^4,  is  a  positive,  real  constant.  R '  is  the  radial 
distance  from  the  reference  point  pf  —  0,  z  =  Zq  to  the  observation  point  ( p\  z  ,  (f)x  )  in  the  far  zone,  and 

the  angle  9X  is  the  elevation  angle  in  the  scatter  direction  counted  from  the  direction  (#0,^.  )  through  the 
observation  point;  see  Figure  5.1. 

Substitution  of  (5.16)  into  the  time-independent  transport  equation  (3.12),  with  the  reduced  incident 
intensity  set  to  zero  since  it  is  very  small  for  Z*  »  1,  pf »  1 ,  gives 


W  rr 

(l-au cos dx )Sca (cos 6X)  =  -2- jj  p(yx )Sn (cos 0[) sin 9[d 0[dy/[ 

4  71  v 


(5.17) 


with 


cos?',  =  cos^,  cos 0X  -I- sin 9X  sin^,'cos(^,  -y/[) 


(5.17a) 


In  Figure  5.1,  only  the  angle  9X  is  shown,  not  the  rotational  azimuthal  angle  y/x  =  <f\  —  (fix  ;  if/x 
could  be  included  in  Figure  5.1  as  the  azimuthal  angle  corresponding  to  the  elevation  angle  9X  defined  in  the 

plane  normal  to  the  direction  9X  =0.  This  definition  for  lj/x  is  allowed  because  the  integral  takes  a  form  which 

is  independent  of  how  the  angular  coordinate  system  is  chosen  since  the  phase  function  depends  only  on  the 
difference  between  the  directions  of  incidence  and  scattering. 

Observe  that  the  integral  equation  (5.17)  for  S^COS#,)  is  identical  to  the  integral  equation  (5.10) 

for  (COS  9  ),  except  for  the  substitution  9  — >  9X ,  which  is  purely  a  mathematical  operation.  Hence,  the 
solution  to  (5. 17)  follows  the  development  for  the  solution  to  (5. 10)  given  above. 

5.4  Composite  Approximate  Solution 

The  first  order  multiple  scattering  solution  of  the  time-independent  transport  equation,  which  is  valid 
at  small  penetration  depths  Z  in  the  forest,  was  shown  to  be  given  by  (5.5)  for  the  collimated  beam  and  by 
(5.7)  for  the  spherical  beam.  At  large  penetration  depths  (z'  »  1)  and  for  fixed,  relatively  small  off-axis 
observation  points,  the  diffuse  intensity  was  shown  to  take  the  asymptotic  form  (5.9)  with  <J^  —  <Jf  min . 


Note  that  both  the  asymptotic  attenuation  rate  and  the  asymptotic  scatter  pattern  S ^  depend  only  on 


31 


Mr)  and  W0,  i.e.,  on  the  parameters  of  the  random  medium,  but  are  largely  independent  of  the  parameters 
of  the  incident  beam.  Using  both  the  first  order  multiple  scattering  solution  and  the  asymptotic  solution,  the 
following  composite  approximate  solution  is  proposed  for  all  Z  and  moderate  values  of  p ' 

ld(p’,z',0,iy) *  rrP{p\z\e^)  =  I? (/>', V*  +  Id° (0, z', 0, 0)S„ (cos 9)[\ - e~fit  ]  (5.18) 
where  /j ° (z\ p\ G,\J/)  is  the  first  order  multiple  scattering  solution  for  the  di Ouse  intensity  given  by  (5.5) 
for  the  collimated  beam  and  by  (5.7)  for  the  spherical  beam.  The  pattern  function  S^COS#)  is  found  as 

discussed  in  Section  (5.3)  and  is  independent  of  the  type  of  incident  wave  involved,  i.e.,  the  pattern  function  is 
the  same  for  the  collimated  and  the  spherical  incident  beam  waves  considered  here.  The  attenuation 

coefficient  /?  is  chosen  to  optimize  agreement  between Id°mp  {p\z\ 9  and  the  rigorous  solutions  for 

I  d  ( p ,  z\9,l//)  given  in  Chapters  3  and  4  for  the  collimated  and  spherical  beams,  respectively. 

Observe  that  according  to  (5.18)  thatld°mp(p\z\6^l//)  «  Id° (p\z\9,l//)  for  z  — >  0  ,  while 

Id°mp (p\z\d.,l//)  &  Id° (0,z\0,0)Sao(cos9)  for  Z  — »oo.  The  latter  limit  yields  estimates  for  the 

parameter  a  and  the  constant  \  since  when  z  — »  oo  and  6  —  0 ,  the  scatter  function  5^(1)  =  1  while 

I d  «  1%  ~  Afi  a™,mnZ  and  Id°mp  «  /Jq( 0,rr,0)  ,  which  is  written  here  to  show  it  is  independent  of  if/  . 

Hence,  the  asymptotic  attenuation  rate  <7^  min  can  be  set  equal  to  the  attenuation  rate  of  the  first  order  solution 

Ij°(  0,  z\ 0)  to  obtain  an  estimate  for  the  parameter  a  while  the  constant  can  be  estimated  by  setting  it 

equal  to  the  magnitude  of  Id°( 0, z\ 0)  using  the  estimated  value  of OC  .  Since  a  rigorous  solution  was  found 

for  the  collimated  beam  (see  Chapter  3),  the  constant  is  determined  more  accurately  by  plotting  the  ratio 

Id  o  /  I™q  versus  z  as  z  —>  oo  for  pr  =  0 ,  0  =  0 ,  where  Id  0  is  the  diffuse  intensity  determined  from 

the  rigorous  solution,  / ^  ^  A)  =  e  ^  s“lce  A/  o  ^  Afo  A)  f°r  z'  — >  oo  . 

5.4.1  Collimated  and  Spherical  Beams 

For  the  collimated  beam,  the  term  Id° i,p\z\09y/)  in  (5.18)  is  given  by  (5.5),  which  in  the 
forward  direction  ( 0  =  0 )  is  independent  of  if/  and  reduces  to 
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(5.19) 


rfV,  2',  0  =  0)  =  erf  P^-Pf-0)e-^lw)2  [e3'  -  e~*  ] 

4  71  w  a 


with  <Jt  =  1  —  dW0  and  p(0)  equal  to  the  phase  function  in  (2.8)  for  y  =  0  =  0  .  For  p*  =  0,  0  =  0  the 
parameter  d  can  be  chosen  using  the  limit  IjOmp(0,z\0,l//)  «  Id°  (0,z\0)  «  for  z  — >  00  .  From 
(5.19)  it  then  follows  that  one  could  choose  T  =  Gt  / min  =  (1  —aJV0)  to  find  an  initial  estimate  ford  . 

For  the  spherical  beam,  the  term  Id° (p\z  ,6,y/)  needed  in  (5.17)  is  given  by  (5.7), 

which  is  then  used  to  find  Id  0  (0,  z\ 0, 0)  .  In  the  forward  direction  (0  =  0)  and  on  the  z  -axis  ( p*  —  0 ), 
it  can  be  shown  that  the  first  order  solution  in  (5.7)  reduces  to  the  following  expression  (which  is  independent 

of  V) 


C(o,Z',o 

'  (4;r)2 


e9*  e* 


(5.20) 


z  +z0  z0 

where  Et  (— x)  is  the  exponential  integral  defined  by 


E,(—x)  =  -  f  e‘—  =  y  +  \nx-  f  (l-e  ')  — ,  (5.20a) 

;  '  ,to  1 

the  Euler-Mascheroni  constant/ =  0.57721 56*  ,  and  D(0)  is  the  directive  gain  in  the  forward 
( 0O  =  0 )  direction. 
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6.  NUMERICAL  RESULTS 


6.1  Introduction 

Numerical  solutions  by  their  very  nature  are  approximate  solutions.  Therefore,  attention  must  be 
paid  to  sources  of  inaccuracies.  In  general,  errors  can  be  minimized  by  increasing  the  trunation  value  for  N  at 
the  expense  of  the  computational  time.  As  in  [2]  and  [13]  unless  stated  otherwise,  the  following  global 

parameters  are  used  in  all  of  the  simulations:  WQ  =  0.75,  A y  ~  0.3,  a  —  0.8,  Tf  =  2,  CCQ  =  4>/5, 

A/ a/  =0.012;  see  [2,13]  for  physical  justification  of  these  assumed  values.  The  following  quantities  vary 

among  different  simulations:  N,  Umax ,  kmax ,  w\  p\  z\  6, if/  . 

6.2  Scatter  or  Phase  Function 

Figure  6.1  shows  the  polar  plot  of  both  the  exact  expression  (2.18)  and  the  truncated  series  expansion 
(3.16)  of  the  phase  function.  This  figure  shows  no  difference  between  the  exact  phase  function  and  the 

approximated  phase  function.  However,  Figure  6.2  shows  the  semilog  plot  of  the  error  6 p  between  the  exact 

and  the  approximated  phase  function,  i.e., 

P exact  P approx 

-  (6.1) 

P exact 

Evidently,  the  error  (6.1)  is  small  when  taking  a  large  number  of  terms  for  the  Legendre  polynomial 
expansion. 

6.3  Convergence  in  Received  Diffuse  Power 

Figures  6.3-6. 7  allows  one  to  observe  the  convergence  in  the  numerical  data  that  is  obtained  for  the 
received  diffuse  power  plotted  versus  time.  By  changing  several  parameters,  particularly  by  increasing  the 
value  of  N,  and  for  different  observation  points,  the  convergence  for  received  deffuse  power  curves  is  shown 
to  improve. 

In  Figure  6.3  for  an  incident  beam  that  is  a  plane  wave,  all  of  the  received  diffuse  power  curves  in 
the  crest  region  overlap  with  one  another  for  different  values  of  N.  However,  in  the  trough  region,  the  received 
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diffuse  power  is  observed  to  be  the  least  accurate  for  N  equals  23  and  the  most  improved  for  N  equal  47. 
Thus,  better  convergence  is  obtained  by  increasing  N. 

Similar  to  the  plane  wave  case,  the  collimated  beam  wave  with  beamwidth  of  ten  in  Figure  6.4 
shows  a  problematic  convergence  in  the  trough  region.  Even  when  N  is  thirty-nine,  the  curve  in  the  trough 
region  does  not  converge  well.  Note,  however,  that  values  of  power  in  the  trough  region  are  very  small.This 
may  contribute  to  the  numerical  inaccuracy  that  is  evident  in  this  region,  which  may  require  using  much 
larger  values  of  N  to  improve  the  series  representation  of  the  solution. 

In  Figure  6.5  for  the  plane  wave  case  ( vv  — >  oo ),  the  received  power  in  the  trough  region  seems  to 
converge  better  for  values  near  0M  =  9°  but  appears  to  require  larger  values  of  N  as  0M  approaches  zero. 
Similar  to  the  plane  wave  case  in  Figure  6.5,  Figure  6.6  shows  that  inaccuracies  appear  in  the  trough  region  as 
0M  approaches  zero.  As  was  shown  in  Figure  6.3,  it  is  expected  that  larger  values  of  N  would  yield  more 
accurate  results  in  the  trough  region.  Since  very  low  power  occurs  in  the  trough  region  and  increasing  N 
necessitates  a  considerable  increase  in  computational  time,  values  of  w  greater  than  10  were  not  considered. 
As  will  be  seen  shortly,  such  inaccuracies  in  the  trough  regions  do  not  occur  for  small  beam  widths,  such  as 
W  =  1 .  Regardless  of  the  different  values  of  N  used,  it  appears  as  shown  in  Figure  6.7  for  the  plane  wave 
case  that  as  the  width  w’  of  the  incident  beam  wave  becomes  smaller,  better  convergence  is  obtained  in  the 
trough  region. 

6.4  Comparisons  with  Quadrature  Method 

An  alternative  solution  to  the  scalar  time-dependent  transport  equation  for  narrow  band,  beam  wave, 
pulse  propagation  in  a  forest  is  presented  in  [13].  In  this  method,  referred  to  as  the  Quadrature  or  Q-method, 
the  Fourier- Bessel  transform  is  used  along  with  the  two-dimensional  Gauss  quadrature  formula  and  an 
eigenvalue-eigenvector  technique,  following  the  procedure  developed  by  Chang  and  Ishimaru  for  CW 
propagation  [7,8].  Comparison  of  results  obtained  using  the  Q-method  and  the  method  presented  here, 
referred  to  as  the  Spherical  Harmonics-Hankel  Transform  or  S -method,  are  presented  in  Figures  6.8  -  6.9.  In 
Figure  6.8,  the  curves  that  are  generated  by  the  Q-  method  (lighter  and  thinner  lines)  lie  very  close  to  the 
curves  that  are  generated  by  the  S-method  (darker  and  thicker  lines).  The  two  set  of  curves  in  Figure  6.9 
display  similar  shapes  and  characteristics.  Values  nearly  match  each  other  over  the  crest  region  and  possess 
the  same  shapes  in  the  trough  region.  The  Q-  method  yields  consistently  lower  values  for  received  power  in 
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the  trough  region  and  in  the  vicinity  of  the  crest  maximum.  Since  the  plane  wave  result  was  obtained  using 
the  method  presented  in  [2],  which  is  derivable  from  the  S-method,  it  is  reasonable  to  assume  that  the  S- 
method  is  more  accurate  than  the  Q-method.  This  was  expected  since  the  Q-method  is  highly  numerical  as 
compared  to  the  S-method.  Since  the  plots  of  the  data  generated  by  two  distinct  methods  essentially  agree, 
credence  is  given  to  results  based  on  the  S-method. 

6.5  Additional  Considerations 

The  truncation  at  Umax  >  12  of  the  series  representation  (2.5a)  for  the  Gaussian  incident  pulses  was 
found  to  yield  values  which  agree  sufficiently  with  the  exact  values  determined  from  (2.2).  In  all  the  graphs, 
umax  was  chosen  to  lie  between  thirteen  and  fifteen  in  order  that  the  received  power  curves  also  be  accurate, 
the  choice  of  umax  being  determined  by  a  minimum  error  criteria.  For  example,  in  Figure  6.10,  the 
normalized  reduced  incident  power  for  different  values  of  Umax  is  plotted,  which  shows  that  Umax  >  12  is 
needed  to  obtain  sufficiently  accurate  results. 

The  k'msx  value  that  is  required  in  the  truncation  of  the  Hankel  transform  in  (3.31)  is  taken  so  that 
more  than  99%  of  the  integrand  is  included.  To  ensure  that  this  is  the  case,  k'max  is  selected  to  be  5  /  w  , 
where  w'  is  the  normalized  beamwidth. 

The  Gaussian  quadrature  method  was  used  to  perform  the  two  integrations  needed  in  the  simulation. 
The  k’-integration  (3.31)  is  approximated  by  thirty-two  terms  while  there  are  ninety-six  terms  for  the 
integration  over  y,  which  is  needed  to  determine  the  expansion  coefficient  g/  in  (3.19). 

The  solutions  to  the  linear  system  of  equations  (3.28)  for  finding  the  particular  solution  and  (3.39) 
for  solving  the  boundary  condition  at  z'  =  0  were  obtained  by  using  LU  factorization,  linear  equation  solver, 
and  iterative  refinement  packages  provided  by  the  optimized  LAPACK  library.  The  eigenvalue  solutions 
were  obtained  by  using  the  ZGGEVX  routine  from  the  LAPACK  library  that  is  based  on  the  QZ  method  in 

EISPACK  [11].  When  tested,  both  of  these  procedures  gave  absolute  errors  ranging  from  10” 10  to  10-16 . 
This  accuracy  is  acceptable,  as  the  library  routines  handle  all  variables  using  double  precision,  which  means  a 
precision  that  is  accurate  up  to  about  sixteen  digits. 

Recall  that  the  phase  function  is  normalized  in  (  2.8a)  such  that 
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(6.2) 


jj  p(y)dQ  =  4  n 

4  n 

This  dictates  that  go  in  (3.19)  must  equal  unity.  However,  go  does  not  equal  unity  when  determined 
numerically  from  (3.19),  which  gives  g0  =  0.9881  when  A ys  =  0.3, a  =  0.8  .  To  ensure  that  g0  *s  unity, 
the  phase  function  p(y)  is  redefined  as  pnorm{/)  =  piy}/  g$  *  which  guarantees  that  (6.2)  is  satisfied  with 
Piy)  replaced  by  pnorm  (y)  . 

6.6  Boundary  Condition 

The  diffuse  intensity  Id  0  is  plotted  versus  0  for  p’  =  0 ,  z'  =  0  and  W  =  1 0  in  Figure  6. 1 1  to 
show  that  the  boundary  condition  (3.2c)  for  the  case  u  =  0  is  satisfied,  i.e.,  that 

7^0=0  at  z  -  0  ,  0<6<nl2  (6.3) 

As  seen  in  the  graph,  when  0  lies  between  0°  and  90°,  the  diffuse  intensity  is  nearly  zero.  As  N  increases,  the 
diffuse  intensity  becomes  even  smaller  and  closer  to  zero  at  z  —  0  over  the  forward  range  0<6<7l/2. 

Figure  6.12  shows  that  for  different  beamwidths  w  ,  the  boundaiy  condition  (6.3)  is  again  very  well 
satisfied  for  the  case  u  =  0  .  Numerical  inaccuracies  produce  the  negligibly  small,  non-zero  values  for  Id  0 

over  the  range  0  <  0  <  90° .  For  u  *  0 ,  the  boundary  condition  Idu=0  at  z'  =  0  ,  0<  0  <  ?r/2  is 
also  well  satisfied  as  shown  in  Figure  6.13  for  o  =  1 .  Observe  that  the  values  of  | Id  ,|  are  significantly 

smaller  than  |/^  0 1  as  indicated  in  Figure  6.13  using  the  same  parameters.  Since  u  *  0 ,  the  diffuse  intensity 
is  complex,  so  its  magnitude  is  plotted  as  distinct  from  the  u  =  0  case,  which  is  purely  real. 

6.7  Power  Attenuation,  Pulse  Broadening  and  Distortion 

In  Figure  6.14,  the  diffuse  intensify  IdQ  at  different  off-axis  locations  initially  increases  rapidly, 

reaches  its  maximum  in  the  vicinity  of  z  -  1  and  then  attenuates  more  slowly  as  the  beam  penetrates  deeper 
into  the  forest.  The  attenuation  is  also  shown  in  the  plots  of  received  power  in  Figures  6.15-6.17  in  the 

forward  direction  (  0M  =  0 )  and  in  Figures  6.18-6.19  in  the  direction  0M  =  4.83°  for  different  values  of 

penetration  depths  z  and  different  off-axis  locations  p' .  The  largest  power  is  transmitted  in  the  forward 
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direction  (0M  =  0)  on  the  beam  axis  (p'  =  0).  When  the  point  of  observation  moves  away  from  the  Z  -axis, 

the  power  decreases.  Distortion  in  the  trough  region  (  see  e.g.  the  curve  in  Figure  6.16  for  p'  =  3  )  occurs 
due  to  small  numerical  inaccuracy  at  the  very  low  power  levels. 

Figure  6.18  clearly  indicates  that  as  the  observation  point  moves  away  from  the  z  =  0  boundary 
between  the  forest  and  the  air,  the  received  power  not  only  attenuates  but  distorts  due  to  pulse  spreading. 
Figure  6. 19  shows  the  received  power  for  different  beamwidths  at  different  off-axis  locations,  but  at  the  same 
penetration  depth.  Observe  that  as  w  increases,  the  received  power  is  stronger,  the  pulse  shapes  remain 
similar,  but  less  distortion  occurs  in  the  trough  region.  The  curves  of  received  power  that  fell  below  -70  dB 
are  not  included  -  because  of  inaccuracies.  Figure  6.20  shows  that  in  the  crest  region,  all  curves  lie  close 
together.  The  only  difference  occurs  in  the  trough  region  where  more  pulse  distortion  occurs  as  the 
beam  width  w  gets  smaller.  Note  that  the  W  =7  case  lies  extremely  close  to  the  plane  wave  result,  which 
shows  the  correct  behavior  of  the  solution  as  the  beamwidth  approaches  large  values. 

Figures  6.21-6.23  depict  power  received  at  larger  scan  angles.  In  Figure  6.21  for  0M  =  62°,  the 
received  power  is  extremely  low  and  exhibits  considerable  distortion  as  compared  to  the  received  power  in 
Figure  6.20  for  0M  =  4.83°  .  The  smaller  widths  possess  the  lowest  received  power  in  the  trough  region. 

The  power  received  for  the  beam  wave  with  w  =7  resembles  the  plane  wave  case  as  noted  previously.  Figure 
6.22  shows  that  the  received  powers  lose  their  distinctive  pulse  shape  for  90 °<0M  <180°.  In  this  range,  the 

time  variation  of  P’  remains  fairly  constant  but  with  magnitudes  that  continue  to  decrease  as  0M  approaches 
180°.  Figure  6.23  shows  what  happens  to  the  received  power  for  the  range  0  <  0M  <  50°  .  Observe  that  the 
crest  of  the  received  power  decreases  as  0M  becomes  larger.  These  graphs  also  indicate  the  effect  of  pulse 

broadening  as  the  beam  penetrates  further  into  the  forest.  Figures  6.21-6.23  indicate  that  as  6 M  increases 
from  0°,  the  antenna  receives  less  power  and  the  pulse  widens,  indicating  stronger  multiscattering,  while  in 
the  scan  directions  90°<6M  <180°,  the  pulse  looses  its  distintive  pulse  shape;  this  happens  because  the 
antenna  now  faces  the  unbounded  region  ( z 9  — >  oo)  in  which  considerable  multi  scattering  occurs. 
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6.8  Angular  Distribution  of  Diffuse  Intensity 

Figure  6.24  and  6.25  depict  the  scatter  patterns  of  diffuse  intensity  Id  0  (time- independent  case)  at 
observation  points  in  the  forest  both  on  and  off  thez'-axis  at  values  of  =0,1, 2,2.5,  and  3  and  at  two 
different  penetration  depths  (z'  =  1  and  5).  In  Figure  6.24  in  the  cross-sectional  plane  defined  by  if/  =0  and 
K  at  the  smaller  value  z  =  1 ,  observe  that  scattering  near  the  axis  is  strongest  in  the  forward  6  =  0 
direction,  but  as  the  observation  point  gets  further  away  from  the  Z  -axis  ( i.e.,  for  larger  values  of  pf ),  side 

scattering  occurs.  For  larger  values  of  z  ,  the  patterns  initially  broaden.  At  the  larger  value  z  =5,  the 
scatter  patterns  exhibit  a  single  major  lobe  that  is  tilted  more  away  from  z  -axis  the  larger  the  value  of  p 9 . 
At  the  observation  point  on  the  z  -axis  ( p'  =  0) ,  the  major  lobe  is  symmetrical  about  its  maximum  in  the 
6  =  0  direction.  Additional  plots  (not  presented  here)  showed  that  at  even  larger  penetration  depths,  at 
observation  points  both  on  and  off  the  z  -axis,  the  scatter  patterns  all  coverged  to  a  symmetrical,  narrow  lobe 
with  a  maximum  in  the  6  =  0  direction,  i.e.  This  observation  led  to  the  development  of  the  asymptotic 
solution  to  the  transport  equation  that  was  presented  in  Chapter  5.  In  Figure  6.25  in  the  cross-sectional  plane 
defined  by  ^  =  —  zr/ 2  and  n/ 2,  the  patterns  initially  broaden,  become  more  isotropic  in  the  backscatter 

direction  in  the  range  Kj2<6  <n ,  and  narrow  in  the  forward  scattering  range  0<6<n/2  (  all  the 
while  diminishing  in  magnitude)  and  eventually  -  though  not  shown  -  coalesces  into  a  single  forward  lobe 
(for  all  values  of  p ' )  with  a  maximum  in  the  6  =  0  direction. 

6.9  Approximate  Solutions 

In  Figures  6.26-6.36,  the  diffuse  intensity  Id  0  for  the  collimated  beam  wave  is  plotted  versus 

normalized  penetration  depth  z  using  the  rigorous  Spherical  Harmonics-Hankel  Transform  (S)  solution 
(3.31)  for  U  =  0,  the  first  order  solution  (5.5a)  and  the  composite  solution  (5.18).  For  p'=0  and  6  =0,  the 
composite  solution  is  identical  to  the  first  order  solution.  For  this  case,  three  solutions  are  plotted  in  Figure 
6.26  for  r  =1,  1.5  and  2,  which  yields  a  =0.716,  0.408  and  0.099,  respectively,  since  T  =  Ot  /  <7^=1- 

dW0  .  From  the  curves,  the  first  order  solution  for  T  =1  yields  the  closest  agreement  with  the  rigorous  S- 
method  solution  over  the  short,  but  important  range  0  <  Z  <  6.5 ,  after  which  the  first  order  solution  for 
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r  =1.5  agrees  best  with  the  rigorous  S  solution  over  the  large  middle  range  6.5  <  z  <  50  .  For  z  >  50,  the 
first  order  solution  for  T  =1  again  gives  the  best  agreement.  Note,  however,  that  the  diffuse  intensity  is  very 
small  for  large  penetration  depths.  The  amplitude  of  the  asymptotic  solution  (5.9),  referred  to  in  the 
figures  as  Asymptotic  #1,  was  determined  as  explained  in  Section  5.4  and  was  used  here  for  both  the 
amplitudes  A0  in  (5.9)  and  Ax  in  (5. 16);  the  latter  is  referred  to  as  Asymptotic  #2  in  the  figures. 

For  the  case  p'  =0  and  9  =  30° ,  Figure  6.27  shows  that  the  composite  solution  for  T  =1.5  agrees 
best  with  the  rigorous  S-method  solution  over  the  range  0  <  z  <  50 ,  after  which  the  composite  solution  for 
r  =1  remains  closest  to  the  S-solution.  For  this  case,  the  values  of  (3  where  chosen  (by  trial  and  error)  such 

that  the  maxima  of  the  composite  solutions  coincide  closely  with  the  maximum  of  the  rigorous  solution;  see 
Figure  6.27a. 

Figure  6.28  for  pf  = 0  and  6  —  30°  shows  that  a  different  choice  for  the  parameter  (3  for  the 
composite  solutions  based  on  optimizing  the  results  in  Figure  6.33,  rather  than  Figure  6.28,  are  not  good 
choices.  Neither  the  first  order  nor  the  composite  solutions  agrees  with  the  S-solution  in  Figure  6.28a. 

Figure  6.29  shows  that  the  composite  solution  for  T  =1  gives  a  good  estimate  for  the  S-solution  for 
0  <  z  <  2  ;  the  first  order  then  is  the  better  one, but  soon  the  composite  solution  for  T  =1.5  is  best  until 
Z  «  30  when  the  first  order  solution  then  again  agrees  better  with  the  S-solution.  In  Figure  6.29  note  that  the 
(3  values  used  were  based  on  optimizing  the  composite  solutions  in  Figure  6.27.  Clearly,  the  composite 
solution  has  limited  utility  when  (3  is  not  optimized.  Similar  comments  apply  to  the  remaining  Figures  6.29- 
6.36. 

From  the  data,  it  was  determined  that  in  the  forward  direction  (0  =  0)  the  first  order  solution 
for  r  =1  can  be  used  to  approximate  the  rigorous  S-solution  and,  therefore,  is  of  use  for  real  time  situations  (it 
does  not  depend  on  the  choice  of  (3 ).  This  is  significant  because  no  information  is  needed  from  the  rigorous 
solution.  For  9  ^  0,  the  first  order  solution  is  again  useful,  but  for  penetration  depths  z’  >  10  .  Small 
penetration  depths  require  adjustment  of  the  attenuation  parameter  (3  to  make  the  composite  solution  useful. 
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For  0^0  and  p  >  0 ,  the  parameter  /?  for  the  composite  solution  needs  to  be  optimized  for  each  value  of 
6  and  p' ,  i.e.,  p  can  be  selected  to  yield  good  results  for  the  composite  solution  for  the  specified  values  of 
6  and  pf ,  but  may  not  give  good  results  for  other  values. 

The  curves  presented  in  Figures  6.26-6.29  serve  to  illustrate  how  approximate  solutions  can  be  used 
to  obtain  sufficiently  accurate  estimates  for  the  intensity  of  the  signal  in  vegetation  in  real  time  for  the  soldier 
in  the  field. 
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7.  CONCLUSION 


The  scalar  time-dependent  equation  of  radiative  transfer  was  used  to  develop  several  theories  of 
pulse  beam  wave  propagation  and  scattering  in  vegetation,  a  medium  characterized  by  many  random  discrete 
scatterers,  which  scatter  energy  strongly  in  the  forward  scattering  direction.  The  specific  problem  analyzed  is 
that  of  a  periodic  sequence  of  Gaussian  shaped  pulses  incident  from  free  space  onto  the  planar  boundary 
surface  of  a  random  medium  half- space,  such  as  a  forest,  that  possesses  a  power  scatter  (phase)  function 
consisting  of  a  strong,  narrow  forward  lobe  superimposed  over  an  isotropic  background.  After  splitting  the 
specific  intensity  into  the  reduced  incident  and  diffuse  intensities,  the  solution  of  the  transport  equation 
expressed  in  cylindrical  coordinates  in  the  random  medium  half-space  was  obtained  by  expanding  the  angular 
dependence  of  both  the  scatter  function  and  the  diffuse  intensity  in  terms  of  Associate  Legendre  polynomials, 
by  using  a  Fourier  series/Hankel  transform  to  obtain  the  equation  of  transfer  for  each  spatial  frequency,  and  by 
satisfying  the  boundary  conditions  that  the  forward  traveling  diffuse  intensity  be  zero  at  the  inter  face  and  zero 
at  infinity. 

Plots  of  intensity  and  received  power  in  the  random  medium  (forest)  showed  distortion  due  to  pulse 
broadening,  angular  spread,  power  attenuation  (especially  at  large  penetration  depths),  and  out-of-beam 
scattering.  Comparison  with  a  second  completely  different  method  of  solution  to  the  scale  transport  equation, 
called  the  Quadrature  method,  substantiated  the  results. 

In  addition,  three  new  approximate  theories  were  developed  for  beam  waves.  The  first  approximate 
theory,  a  first  order  multiple  scattering  theory,  allowed  for  the  estimation  of  the  intensity  at  small  penetration 
depths  without  resorting  to  computer  intensive  calculations  while  maintaining  good  accuracy.  The  second 
new  theory  yielded  an  asymptotic  solution  to  the  transport  equation  which  is  valid  at  large  penetration  depths 
into  the  forest.  This  allowed  a  characterization  of  the  field  deep  in  the  forest  also  without  extensive 
calculations.  The  third  approximate  theory  was  a  composite  solution  which  combined  both  the  first  order 
solution  and  the  asymptotic  solution  and  closed  the  gap  between  the  first  two  approximate  theories. 

The  three  approximate  theories  were  numerically  compared  to  the  rigorous  theory.  The  first  order 
theory  was  shown  to  agree  with  the  rigorous  theory  at  short  vegetation  penetration  depths.  At  larger  depths,  it 
agrees  in  the  forward  scatter  direction  only,  but  not  otherwise.  The  asymptotic  theory  was  shown  to  have  the 
correct  behavior  in  all  scatter  directions  and  to  agree  with  the  rigorous  theory  at  large  penetration  depths.  The 
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development  of  an  asymptotic  theory  is  new;  its  feasibility  was  predicted  by  the  extensive  numerical  data  that 
was  obtained  from  the  rigorous  theory.  One  important  aspect  of  the  asymptotic  theory  is  that  it  is  independent 
of  the  incident  radiation  and,  therefore,  remains  valid  for  other  cases  as  well  and,  in  particular,  for  the 
spherical  beam  case. 

As  opposed  to  the  rigorous  theory,  the  three  approximate  theories  are  numerically  efficient.  The 
significance  of  this  is  that  they  are  useful  for  near-real-time  modeling  of  the  wireless  microwave 
communication  channel  in  vegetation  and,  therefore,  of  interest  to  the  soldier  in  the  field. 
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Figure  2.1  Bounded  beam  wave  pulse  train  incident  onto  a  forest  half-space 

(a)  Collimated  incident  beam  wave  pulse  train  with  p  -dependent  Gaussian  amplitude  A  (p) 

(b)  Spherically  diverging  incident  beam  wave  pulse  train  with  directive  gain  Z)  (#0  )  and  the 
phase  center  of  the  beam  located  at  Z  =  ZQ  <  0 
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Figure  2.2  The  basic  geometry  for  scattering  in  the  forest  half-space  of  an  incident  beam  wave  pulse 
train.  The  “sphere”  represents  a  scatter  point  in  the  forest.  The  tilted  cylinder  represents  a 
receiving  antenna,  which  is  shown  enlarged  in  Figure  2.3. 
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Figure  23 


The  coordinate  geometry  of  the  receiving  antenna,  which  is  depicted  as  a  tilted  cylinder  with  a 
main  beam  direction  (dM ,  y/ M  ) ;  note  0  =  0M  . 
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Figure  6.2 


Error  analysis  of  the  phase  function,  where  N’=N+1  is  the  number  of  terms  for  the  truncation 
of  (3. 16). 
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Figure  63 


Normalized  received  diffuse  power  Pd 
for  w  — >  oo  (plane  wave),  p  =0, 


o  =15. 

max 


versus  normalized  time  t'  for  different  values  of  N 
z'  =  1 ,  Gy  =0  ( y/M  is  undefined  when  pf  =  0  ), 
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Figure  6.4  Normalized  received  diffuse  power  Pd’  versus  normalized  time  t'  for  N=23,  27,  39  and 
w'=10,p'  =  0  ,z'  =  \,0M  =  0 ,  vmax  =15. 
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Figure  6.5 


Normalized  received  diffuse  power  Pd  versus  normalized  time/' for  different  values  of  0M 
in  the  range  0°  <  6M  <  9°  and  w'  — >  oo  (plane  wave),  p'  =  0 ,  z  =  1 ,  L>max  =15,  N=3 1 . 
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Normalized  Time,  t’ 


Figure  6.6  Normalized  received  diffuse  power  Pd’  versus  normalized  time/' for  w  =  10,  p’=0,  z’=l, 
Gm  =0°  and  4.83°,  Umax  =15,  N=39. 
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Figure  6.7  Normalized  received  diffuse  power  Pd*  versus  normalized  timef'  for  w'  =  1,10 
andw'  — »  oo  (plane  wave)  with  associated  values  of  N  =  31,  39  and  47,  respectively,  and 
P’=0>  z’=l,  0M  =0,t>max  =  15 . 
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Figure  6.8  Comparisons  between  Quadrature  method  (Q)  and  Spherical  Harmonics- Hankel  Transform 
method  (S)  of  normalized  received  power  P’  versus  normalized  time  t’  for  z’=l,3,5,10  and 

w'  =  1 ,  p’=0,  eM  =4.83 °,L>max  =  1 5 ,  N=3 1 . 
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Figure  6.9  Comparison  between  Quadrature  method  (Q)  and  Spherical  Harmonics-Hankel  Transform 
method  (S)  of  the  normalized  received  power  P’  versus  normalized  time  t’  for 

w'  =  0.5,l,3,5,oo  and  p’=0,  z’=l,  0M  =4.83°,L>max  =15,  N=27. 
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Normalized  Received  Reduced  Incident  Power,  P‘ri[dB] 


Figure  6.10  Normalized  reduced  incident  power  P\  versus  normalized  time  t’  for  p’=0,  z’=l, 
different  values  of t>max . 
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Diffuse  Intensity,  Id,0 


Figure  6.11  Diffuse  Intensity  Id  0  versus  0  for  w  =10,  p'  =  0,z'  =  0 ,  u=0  and  N=23,27,31. 
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Figure  6.12  Diffuse  intensity  Id0  versus  6  for  beam  widths  w'  =  0.5,1,2,3,5,7,10,  cc  and  p’=0,  z’=0, 
u=0,  N=31. 
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Magnitude  of  diffuse  intensity  \Id  ,  versus  6  for  w  =1,  2  and  p’=0,  z’=0.  o=l,  N=31. 
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Figure  6.14 


Diffuse  Intensity  I d  0  versus  normalized  penetration  depth  z’  for  \v  =  1 ,  p’=0,  0=0,  u=0. 


N=31. 
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Figure  6.15  Normalized  received  power  P’  versus  normalized  time  t’  for  p  varying  from 
and  w'  =  \,  z’=l,  GM  =0,  [// M  =0,  t>max=15,  N=31. 
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Normalized  Received  Power,  P’  [dB] 


Figure  6.16  Normalized  received  power  P’  versus  normalized  time  t’  for  p '  =  0,1, 2, 3  and  W  =  1 ,  z’=3, 
Vm  =0,umax=15,N=27. 
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Normalized  Received  Pcwer,  P'  icB]  Normalized  Received  Power,  p1  icbi 

i  i  i  i  i  i  i 


Figure  6.17  Normalized  received  power  P’  versus  normalized  time  t’  for  z’=0.5,0.7,l,1.5,2,3 
and  W  =  1 , 0U  =0,^  =0,  umax  =15,  N=27,  (a)  p'  =0 ,  (b)  p'  =  1 ,  (c)  p'  =  2,  (d)  p' =  3  . 
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Figure  6.18  Normalized  received  power  P’  versus  normalized  time  t’  for  z’=l,3,5,10  and  w'  =  1 ,  p’=0, 
^=4.83°,  y/M  =0,  umax=  15.N-27. 
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Figure  6.19  Normalized  received  power  P1  versus  normalized  time  t’  for  (a)  vt/=2,  (b)  w/ =3,  (c)  w'  =5, 
(d)  w'=7  with  p’=0,  1,5  ,  10  and  for  z’=l,  0M  =4.83°,  y/M  =0,  ^ax  =15,  N=27. 
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Figure  6.20  Normalized  received  power  versus  normalized  time  t’  for  collimated  beam  wave  with 
beamwidths  =0.5, 1,2, 3, 5, 7  and  for  a  plane  wave  (wf— >oo),  p=0,  z’=l, 

0M  =4.83°,  ^max  =15;  N=31  for  the  plane  wave  case  and  N=27  for  collimated  beam  waves. 
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Figure  6.21  Normalized  received  power  P'  versus  normalized  time  t’  for  collimated  beam  waves  with 
beamwidths w'=0.5,  1,  2,  3,  5,  7  and  for  a  plane  wave  (w— >oo),  p=0,  z’=l, 
0M  =62°,  L>max  =  15,  N=31  for  the  plane  wave  case  and  N=27  for  collimated  beam  waves. 
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Figure  6.22  Normalized  received  power  P'  versus  normalized  time  t’  for  collimated  beam  waves  with 
beamwidths  vv  =0.5,1, 2,3,5,7  and  for  a  plane  wave  ( W  — >  oo ),  p=0,  z’=l,  (a)  dM  =87°,  (b) 

0M  =1 18°,  (c)  9m  =150°,  (d)  6m  =175°,  t>max  =  15,  N=31  for  the  plane  wave  case  and  N=27 
for  collimated  beam  waves. 
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Normalized  Received  Pcwer,  P’  icBi  Normalized  Received  Pa/ver,  P'  icB] 
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Figure  6.23  Normalized  received  power  P’  versus  normalized  time  t’  for  different  values  of  dM  and  for 
(a)  z’=l,  (b)  z’=3,  (c)  z’=5,  and  (d)  z’=10  with  wf= 1,  p’=0,  umax  =  15,  N=27. 
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Figure  6.24  Diffuse  Intensity  Id0  vs  6  for  p'  =  0,1, 2, 2.5, 3  at  (a)  z’=l  and  (b)  z’=5  with  vv' 
if/  =  0,71  andN=31. 
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Figure  6.25  Diffuse  Intensity  Id0  vs  0  for  p'  =  0,1, 2, 2.5, 3  at  (a)  z’=l  and  (b 
\f/  =-nj2,nl2  andN=31. 
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Figure  6.26  Diffuse  Intensity  Id  0  vs  normalized  penetration  depth  Z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p'  =  0  and  6  =  0°  for  (a)  0  <  z  <  10  and  (b) 
0  <  z  <  100  with  w=  1,  [//  =  0  and  N=31.  Note  that  the  values  of  /?  do  not  apply  in  this 
case. 


74 


(a) 


(b) 

Figure  6.27  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p9  =  0  and  6  =  30°  for  (a)  0  <  Z  <  1 0  and  (b) 
0<z'<100  with  vi/  =  l,  l//  =  0  and  N=31.  Note  that  the  values  of  /?  are  chosen  to 
optimize  the  composite  solution  in  this  case. 
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Figure  6.28  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p'  =  0  and  0  =  30°  for  (a)  0  <  Z  <  1 0  and  (b) 
0<zf<100  withw,  =  l,  l// =  0  and  N=31.  Note  that  the  values  of  /?  are  chosen  to 
optimize  the  composite  solution  for  p  —  2  and  6  =  0°  in  Figure  6.34. 
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Figure  6.29  Diffuse  Intensity  Id  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p'  =  1  and  0  =  0°  for  (a)  0  <  z  <  10  and  (b) 
0<z'<100  with  W =\,  l//  =  0  and  N=31.  Note  that  the  values  of  (3  are  chosen  to 
optimize  the  composite  solution  for  p'  =  0  and  6  =  30°  in  Figure  6.27. 
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Figure  630  Diffuse  Intensity  1  d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  pf  =  1  and  9  =  0°  for  (a)  0  <  z'  <  10  and  (b) 
0<z'<100  withw=l,  l//  =  0  and  N=31.  Note  that  the  values  of  /?  are  chosen  to 
optimize  the  composite  solution  for  p'  =  2  and  6  —  0°  in  Figure  6.34. 
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Figure  631  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p’  =  1  and  0  =  30°  for  (a)  0  <  z  <  1 0  and  (b) 
0<z'<100  withw/=l,  if/ =  0  and  N=31.  Note  that  the  values  of  [5  are  chosen  to 
optimize  the  composite  solution  for  p'  =  0  and  6  =  30°  in  Figure  6.27. 
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Figure  632  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p  =  1  and  6  —  30°  for  (a)  0  <  z  <  1 0  and  (b) 
0<z'<100  withw,=l,  l// =  0  and  N=31.  Note  that  the  values  of  /?  are  chosen  to 
optimize  the  composite  solution  for  p'  =  2  and  6  =  0°  in  Figure  6.34. 
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Figure  6 33  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p'  =  2  and  9  —  0°  for  (a)  0  <  Z  <  10  and  (b) 
0<z'<100  with  w  =1,  l//  =  0  and  N=31.  Note  that  the  values  of  (5  are  chosen  to 
optimize  the  composite  solution  for  p'  =  0  and  6  =  30°  in  Figure  6.27. 
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Figure  634  Diffuse  Intensity  Id  0  vs  normalized  penetration  depth  Z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p9  =  2  and  6  —  0°  for  (a)  0  <  z  <  10  and  (b) 
0<zf<100  withw'  =  l,  l//  =  0  and  N=31.  Note  that  the  values  of  (5  are  chosen  to 
optimize  the  composite  solution  in  this  case. 
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Figure  635  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  pf  =  2  and  6  =  30°  for  (a)  0  <  z  <  10  and  (b) 
0<z'<100  withwf=l,  l//  =  0  and  N=31.  Note  that  the  values  of  (3  are  chosen  to 
optimize  the  composite  solution  for  p'  =  0  and  6  =  30°  in  Figure  6.27. 
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Figure  6 36  Diffuse  Intensity  I d  0  vs  normalized  penetration  depth  z  for  the  collimated  beam  case  using 
the  rigorous  and  approximate  solutions  for  p'  =  2  and  6  =  30°  for  (a)  0  <  z  <  10  and  (b) 
0<z'<100  withwf=l,  l//  =  0  and  N=31.  Note  that  the  values  of  /?  are  chosen  to 
optimize  the  composite  solution  for  p'  =  2  and  6  =  0°  in  Figure  6.34. 
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